Práctica Final
Alejandro Cárdenas Barranco, Álvaro Rodríguez Gallardo, Juan Manuel Rodríguez Gómez
2024-01-05
En esta práctica se va a analizar un dataset que aborda el rendimiento en la asignatura Matemáticas de los alumnos en educación secundaria de dos centros escolares portugueses. Los datos se han recopilado mediante el uso de informes y cuestionarios escolares. Dicho dataset se puede descargar mediante el siguiente enlace:
https://archive.ics.uci.edu/dataset/320/student+performance
En primer lugar, se realizará un análisis exploratorio previo de los datos para identificar posibles valores perdidos y valores extremos o “outliers”, y se tomarán las decisiones correspondientes para tratarlos.
En segundo lugar, se realizará un Análisis de Componentes Principales (ACP), esto es, se va a condensar la información aportada por las variables originales en unas pocas combinaciones lineales de ellas, con el objetivo de hacer una reducción de la dimensión. Estas se buscan con máxima varianza y perpendiculares entre sí, de forma que cada una sigue la dirección en la que las observaciones varían más y no está correlacionada con las anteriores.
En tercer lugar, se realizará un Análisis Factorial (AF), esto es, se van a identificar variables latentes (no observables) con una alta correlación con un grupo de variables observables y correlación prácticamente nula con el resto. Así, se hará una reducción de la dimensión.
En cuarto lugar, se realizará un Análisis Discriminante Lineal (ADL) y otro Cuadrático (ADC), verificando las hipótesis necesarias de normalidad. El Análisis Discriminante es un método de clasificación de variables cualitativas que permitirá clasificar nuevas observaciones según sus características (variables explicativas o predictores) en las distintas categorías de la variable cualitativa respuesta.
Finalmente, se realizará un Análisis Cluster (AC). El Análisis Cluster es una técnica multivariante cuyo principal objetivo es agrupar objetos formando conglomerados (clusters) con un alto grado de homogeneidad interna y heterogeneidad externa. En otras palabras, el AC es un procedimiento exploratorio que permite encontrar estructuras con similitudes en un conjunto de datos con cierta variabilidad.
1 Cargando Paquetes y Leyendo el Conjunto de Datos
1.1 Cargar e instalar paquetes de R
El siguiente módulo de código fuente se encarga de cargar, si ya están instalados, todos los paquetes que se utilizarán en esta sesión de R. Si bien un paquete R se puede cargar en cualquier momento cuando se vaya a utilizar, es recomendable optimizar sus llamadas con este fragmento de código al principio.
Cargar un paquete en una sesión de R requiere que ya esté instalado. Si no es así, el primer paso es ejecutar la sentencia:
install.packages(“name_of_the_library”)
#########################################
# Loading necessary packages and reason #
#########################################
# This is an example of the first installation of a package
# Only runs once if the package is not installed
# Once it is installed this sentence has to be commented (not to run again)
# install.packages("summarytools")
# Package required to call 'read_xlsx' function (loading '.xlsx' data format)
library(readxl)
# Package required to call 'freq' and 'descr' functions (descriptive statistics)
library(summarytools)
# Package required to call 'factanal' function
library(stats)
# Package required to call 'stat.desc' function
library(pastecs)
# Package required to call 'cortest.bartlett' function
library(psych)
# Package required to call 'hetcor' function
library(polycor)
# Package required to call 'ggcorrplot' function
library(ggcorrplot)
# Package required to call 'corrplot' function
library(corrplot)
# Package required to call 'rplot' function
library(corrr)
# Package required to call 'ggplot' function (graphical tools)
library(ggplot2)
# Package required to call 'ggarrange' function (graphical tools)
library(ggpubr)
# Package required to call 'fviz_pca_var, fviz_pca_ind and fviz_pca' functions
library(factoextra)
# Package required to call 'scatterplot3d' function
library(scatterplot3d)
# Package required to call 'mutate' function
library(tidyverse)
# Package required to call 'clusGap' function
library(cluster)
# Package required to call 'ggdendrogram' function
library(ggdendro)
# Package required to call 'grid.arrange' function
library(gridExtra)
# Package required to call 'melt' function
library(reshape2)
# Package required to call 'mvn' function
library(MVN)
# Package required to call 'boxM' function
library(biotools)
# Package required to call 'partimat' function
library(klaR)
# Package required to call 'summarise' function
library(dplyr)
# Package required to call 'createDataPartition' function
library(caret)1.2 Descripción del conjunto de datos student_mat_DB
El fichero student_mat_DB contiene datos recogidos de 395 estudiantes recogidos en 31 variables. A continuación se muestra una descripción de cada variable del conjunto de datos:
- Centro escolar del estudiante (school)
- Sexo del estudiante (sex)
- Edad del estudiante (age)
- Tipo de dirección de vivienda del estudiante (address)
- Número de miembros de la familia del estudiante (famsize)
- Estado de convivencia de los padres del estudiante (Pstatus)
- Educación de la madre del estudiante (Medu)
- Educación del padre del estudiante (Fedu)
- Trabajo de la madre del estudiante (Mjob)
- Trabajo del padre del estudiante (Fjob)
- Motivo del estudiante por el cual eligió su centro escolar (reason)
- Tutor del estudiante (guardian)
- Tiempo de viaje de casa del estudiante al centro escolar (traveltime)
- Tiempo de estudio semanal del estudiante (studytime)
- Número de materias suspendidas en el pasado (failures)
- El estudiante recibe apoyo educativo adicional (schoolsup)
- El estudiante recibe apoyo educativo familiar (famsup)
- El estudiante recibe clases extras pagadas de la asignatura Matemáticas (paid)
- El estudiante está apuntado en actividades extracurriculares (activities)
- El estudiante asistió de pequeño a la guardería (nursery)
- El estudiante quiere cursar estudios superiores (higher)
- El estudiante tiene acceso a Internet en casa (internet)
- El estudiante se encuentra en una relación de pareja (romantic)
- Calidad de las relaciones familiares del estudiante (famrel)
- Cantidad de tiempo libre que posee el estudiante después del colegio (freetime)
- Frecuencia con la que el estudiante sale a la calle con sus amigos (goout)
- Cantidad de consumo de alcohol del estudiante en jornada laboral (Dalc)
- Cantidad de consumo de alcohol del estudiante en el fin de semana (Walc)
- Estado de salud actual del estudiante (health)
- Número de veces que el estudiante ha faltado al colegio (absences)
- Nota final del estudiante en la asignatura Matemáticas (G3)
1.2.1 Lectura del conjunto de datos
El archivo de datos tiene extensión .xlsx (Microsoft Excel). Para cargar el fichero se utiliza la función read_excel() dentro del paquete readxl.
Para cargar la base de datos, la sesión de R debe ejecutarse en el mismo directorio en el que se encuentra el fichero. Para ello, se puede usar la función setwd(“Ubicación del archivo de datos”) o se puede pulsar Session/Set Working Directory/Choose Directory en la barra de menús de RStudio. Asumiendo que la sesión de R está correctamente direccionada, el siguiente código carga el fichero de datos en un data.frame.
# Setting the work directory
setwd("C:/Users/LENOVO/Desktop/EMV/practicas/practica_final")
# Loading a .xlsx (excel) file.
# The output of this function is already a data.frame object
# Remember that package 'readxl' is required
data_xlsx<-read_excel("student_mat_DB.xlsx", sheet = 2)
# This sentence identifies the type of object that the identifier represents
class(data_xlsx)## [1] "tbl_df" "tbl" "data.frame"2 Análisis Explotario Univariante
2.1 Recodificaciones o agrupaciones de datos
Se ha recodificado el conjunto de datos original. Esta recodificación ya se encuentra realizada y detallada en el archivo student_mat_DB.xlsx.
Cabe destacar que es cierto que, por el carácter de las variables, principalmente categóricas pero recodificadas, podrían discutirse los resultados obtenidos tras aplicar las técnicas para la reducción de la dimensión (ACP y AF, los cuales realizaremos posteriormente) puesto que estos trabajan con variables cuantitativas, pero también es cierto que si las variables son binarias (como ocurre en varias variables de nuestro conjunto de datos), todavía se puede justificar su uso y usabilidad.
Hay que tener cuidado con las variables codificadas que eran categóricas originalmente y que tienen más de dos niveles de respuesta porque quizá no interese incluirlas. Para decidir si las incluímos o no, tendría que ser el que ha diseñado el dataset el que nos dijera si esas variables son relevantes para el problema bajo estudio. Asumiendo que lo son, una forma de ver si pueden ser útiles es hacer un análisis descriptivo de ellas y ver si los niveles están compensados (más o menos hay el mismo número de individuos en cada uno). Si lo están es una variable que podría aportar información, pero si hay algún nivel muy bajo o si todo se acumula en un nivel, habría que plantearse descartar la variable. No hay una forma definitiva de actuar, depende de lo que busquemos, pero este camino podría ayudar.
2.1.1 Variable Medu
# Frequency tables. Descriptive analysis
# Remember that package 'summarytools' is required
freq(as.factor(data_xlsx$Medu))## Frequencies
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
## 0 3 0.76 0.76 0.76 0.76
## 1 59 14.94 15.70 14.94 15.70
## 2 103 26.08 41.77 26.08 41.77
## 3 99 25.06 66.84 25.06 66.84
## 4 131 33.16 100.00 33.16 100.00
## <NA> 0 0.00 100.00
## Total 395 100.00 100.00 100.00 100.00# Pie chart and bar graph
p1<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$Medu)))+geom_bar()+
coord_polar("y")+labs(x="Medu", y="%")
p2<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$Medu)))+geom_bar()+
labs(x="Medu", y="%")
# This function controls the graphical output device
# Remember that package 'ggpubr' is required
ggarrange(p1, p2, nrow = 1, ncol=2, common.legend = TRUE)
Vemos que la frecuencia del nivel asociado al valor 0 es muy baja (no están compensados los niveles). Por tanto, eliminamos dicha variable de nuestro conjunto de datos.
# The variable 'Medu' is eliminated
data_xlsx<-data_xlsx[, -c(7)]2.1.2 Variable Fedu
# Frequency tables. Descriptive analysis
# Remember that package 'summarytools' is required
freq(as.factor(data_xlsx$Fedu))## Frequencies
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
## 0 1 0.26 0.26 0.25 0.25
## 1 77 20.26 20.53 19.49 19.75
## 2 114 30.00 50.53 28.86 48.61
## 3 97 25.53 76.05 24.56 73.16
## 4 91 23.95 100.00 23.04 96.20
## <NA> 15 3.80 100.00
## Total 395 100.00 100.00 100.00 100.00# Pie chart and bar graph
p3<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$Fedu)))+geom_bar()+
coord_polar("y")+labs(x="Fedu", y="%")
p4<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$Fedu)))+geom_bar()+
labs(x="Fedu", y="%")
# This function controls the graphical output device
# Remember that package 'ggpubr' is required
ggarrange(p3, p4, nrow = 1, ncol=2, common.legend = TRUE)
Vemos que la frecuencia del nivel asociado al valor 0 es muy baja (no están compensados los niveles). Por tanto, eliminamos dicha variable de nuestro conjunto de datos.
# The variable 'Fedu' is eliminated
data_xlsx<-data_xlsx[, -c(7)]2.1.3 Variable Mjob
# Frequency tables. Descriptive analysis
# Remember that package 'summarytools' is required
freq(as.factor(data_xlsx$Mjob))## Frequencies
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
## 0 55 14.29 14.29 13.92 13.92
## 1 33 8.57 22.86 8.35 22.28
## 2 100 25.97 48.83 25.32 47.59
## 3 59 15.32 64.16 14.94 62.53
## 4 138 35.84 100.00 34.94 97.47
## <NA> 10 2.53 100.00
## Total 395 100.00 100.00 100.00 100.00# Pie chart and bar graph
p5<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$Mjob)))+geom_bar()+
coord_polar("y")+labs(x="Mjob", y="%")
p6<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$Mjob)))+geom_bar()+
labs(x="Mjob", y="%")
# This function controls the graphical output device
# Remember that package 'ggpubr' is required
ggarrange(p5, p6, nrow = 1, ncol=2, common.legend = TRUE)
Vemos que más o menos están compensados los niveles de esta variable. Por tanto, mantenemos dicha variable en nuestro conjunto de datos.
2.1.4 Variable Fjob
# Frequency tables. Descriptive analysis
# Remember that package 'summarytools' is required
freq(as.factor(data_xlsx$Fjob))## Frequencies
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
## 0 29 7.34 7.34 7.34 7.34
## 1 18 4.56 11.90 4.56 11.90
## 2 111 28.10 40.00 28.10 40.00
## 3 20 5.06 45.06 5.06 45.06
## 4 217 54.94 100.00 54.94 100.00
## <NA> 0 0.00 100.00
## Total 395 100.00 100.00 100.00 100.00# Pie chart and bar graph
p7<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$Fjob)))+geom_bar()+
coord_polar("y")+labs(x="Fjob", y="%")
p8<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$Fjob)))+geom_bar()+
labs(x="Fjob", y="%")
# This function controls the graphical output device
# Remember that package 'ggpubr' is required
ggarrange(p7, p8, nrow = 1, ncol=2, common.legend = TRUE)
Vemos que la frecuencia del nivel asociado al valor 4 es muy alta (no están compensados los niveles). Por tanto, eliminamos dicha variable de nuestro conjunto de datos.
# The variable 'Fjob' is eliminated
data_xlsx<-data_xlsx[, -c(8)]2.1.5 Variable reason
# Frequency tables. Descriptive analysis
# Remember that package 'summarytools' is required
freq(as.factor(data_xlsx$reason))## Frequencies
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
## 0 109 27.59 27.59 27.59 27.59
## 1 105 26.58 54.18 26.58 54.18
## 2 145 36.71 90.89 36.71 90.89
## 3 36 9.11 100.00 9.11 100.00
## <NA> 0 0.00 100.00
## Total 395 100.00 100.00 100.00 100.00# Pie chart and bar graph
p9<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$reason)))+geom_bar()+
coord_polar("y")+labs(x="reason", y="%")
p10<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$reason)))+geom_bar()+
labs(x="reason", y="%")
# This function controls the graphical output device
# Remember that package 'ggpubr' is required
ggarrange(p9, p10, nrow = 1, ncol=2, common.legend = TRUE)
Vemos que más o menos están compensados los niveles de esta variable. Por tanto, mantenemos dicha variable en nuestro conjunto de datos.
2.1.6 Variable guardian
# Frequency tables. Descriptive analysis
# Remember that package 'summarytools' is required
freq(as.factor(data_xlsx$guardian))## Frequencies
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
## 0 273 69.11 69.11 69.11 69.11
## 1 90 22.78 91.90 22.78 91.90
## 2 32 8.10 100.00 8.10 100.00
## <NA> 0 0.00 100.00
## Total 395 100.00 100.00 100.00 100.00# Pie chart and bar graph
p11<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$guardian)))+geom_bar()+
coord_polar("y")+labs(x="guardian", y="%")
p12<-ggplot(data_xlsx, aes(x=factor(1),fill=as.factor(data_xlsx$guardian)))+geom_bar()+
labs(x="guardian", y="%")
# This function controls the graphical output device
# Remember that package 'ggpubr' is required
ggarrange(p11, p12, nrow = 1, ncol=2, common.legend = TRUE)
Vemos que la frecuencia del nivel asociado al valor 0 es muy alta (no están compensados los niveles). Por tanto, eliminamos dicha variable de nuestro conjunto de datos.
# The variable 'guardian' is eliminated
data_xlsx<-data_xlsx[, -c(9)]2.2 Explorando el conjunto de datos
# Retrieving the name of all variables
colnames(data_xlsx)## [1] "school" "sex" "age" "address" "famsize"
## [6] "Pstatus" "Mjob" "reason" "traveltime" "studytime"
## [11] "failures" "schoolsup" "famsup" "paid" "activities"
## [16] "nursery" "higher" "internet" "romantic" "famrel"
## [21] "freetime" "goout" "Dalc" "Walc" "health"
## [26] "absences" "G3"# Displaying a few records
head(data_xlsx, n=10)data_xlsx# Checking the type of the variables
sapply(data_xlsx, class)## school sex age address famsize Pstatus Mjob
## "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## reason traveltime studytime failures schoolsup famsup paid
## "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## activities nursery higher internet romantic famrel freetime
## "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## goout Dalc Walc health absences G3
## "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"2.3 Valores perdidos (NA)
# Construction of the function that calculates the percentage of missing values for each variable
porcentaje_NA<-function(data){
c=(sum(is.na(data)))/length(data)*100
return(c)
}
# Checking for missing values
cont_NA<-apply(data_xlsx, 2, porcentaje_NA)
cont_NA## school sex age address famsize Pstatus Mjob
## 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.531646
## reason traveltime studytime failures schoolsup famsup paid
## 0.000000 0.000000 13.670886 0.000000 0.000000 0.000000 0.000000
## activities nursery higher internet romantic famrel freetime
## 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## goout Dalc Walc health absences G3
## 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000Se observa que la variable Mjob tiene menos de un 5% de valores perdidos, a diferencia de la variable studytime, que sí que tiene más de un 5%. Por lo tanto, hay que analizar en dicha variable studytime si el patrón seguido por dichos valores perdidos es aleatorio o no. Para ello, se estudia la homogeneidad según grupos (NA y no NA) con otras variables según un Test de Student (ya que la variable es continua).
Se define una nueva variable studytime_NA que toma el valor 0 si el correspondiente dato no es perdido y el valor 1 si el correspondiente dato es perdido. A continuación, se realiza un test de contraste de igualdad de medias entre los grupos de la variable traveltime (que no tiene valores perdidos) correspondientes al 0 y al 1 de estas nuevas variables.
# Construction of the function that replaces missing values of a variable with 1 and non-missing values with 0
indicadora_NA<-function(data, na.rm=F){
data[!is.na(data)]<-0
data[is.na(data)]<-1
return(data)
}
studytime_NA<-indicadora_NA(data_xlsx$studytime)
# Test de Student
t.test(data_xlsx$traveltime~studytime_NA, var.equal=TRUE)##
## Two Sample t-test
##
## data: data_xlsx$traveltime by studytime_NA
## t = -1.2192, df = 393, p-value = 0.2235
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## -0.32519085 0.07624983
## sample estimates:
## mean in group 0 mean in group 1
## 1.431085 1.555556Como el p-valor es grande en todos los casos (p-valor > 0.15), no se encuentran evidencias para rechazar la hipótesis nula de igualdad de medias, de forma que se acepta la hipótesis de homogeneidad y se concluye que el patrón es aleatorio. En este caso, como las variables son cuantitativas, la decisión para los datos perdidos es reemplazarlos por la media de su variable.
# Construction of the function that handles missing values
not_available<-function(data, na.rm=F){
data[is.na(data)]<-mean(data, na.rm=T)
return(data)
}
data_xlsx$Mjob<-not_available(data_xlsx$Mjob)
data_xlsx$studytime<-not_available(data_xlsx$studytime)
# Viewing the data again
head(data_xlsx, n=10)# Checking for missing values
no_cont_NA<-apply(data_xlsx, 2, porcentaje_NA)
no_cont_NA## school sex age address famsize Pstatus Mjob
## 0 0 0 0 0 0 0
## reason traveltime studytime failures schoolsup famsup paid
## 0 0 0 0 0 0 0
## activities nursery higher internet romantic famrel freetime
## 0 0 0 0 0 0 0
## goout Dalc Walc health absences G3
## 0 0 0 0 0 02.4 Análisis descriptivo numérico clásico
A continuación, se realiza un análisis descriptivo de los datos. La función stat.desc() del paquete pastecs proporciona las medidas de posición, dispersión y forma más importantes:
- Número de valores (nbr.val)
- Número de valores nulos (nbr.null)
- Número de valores perdidos (nbr.na)
- Valor mínimo de la variable (min)
- Valor máximo de la variable (max)
- Rango de valores de la variables (max-min) (range)
- Suma de los valores no perdidos (sum).
- Mediana (median).
- Media (mean).
- Error estándar de la media (SE.mean)
- Intervalo de confianza de la media (CI.mean)
- Varianza (var)
- Desviación típiaca (std.dev)
- Coeficiente de variación (coef.var)
- Coeficiente de asimetría (skewness)
- Estadístico para contrastar si el coeficiente de asimetría es cero (skew.2SE)
- Coeficiente de curtosis (kurtosis)
- Estadístico para contrastar si el coeficiente de curtosis es cero (kurt.2SE).
- Estadístico de Shapiro-Wilks (normtest.W)
- Probabilidad asociada al estadístico de Shapiro-Wilks (normtest.p)
# Basic descriptive statistics
stat.desc(data_xlsx, norm=TRUE)Para más detalles sobre esta función se puede consultar el siguiente enlace:
https://www.rdocumentation.org/packages/pastecs/versions/1.3.21/topics/stat.desc
Como esta función no muestra los cuartiles, se utiliza la función summary() para obtenerlos junto al resto de medidas de posición:
# Basic descriptive statistics
summary(data_xlsx)## school sex age address
## Min. :0.0000 Min. :0.0000 Min. :15.0 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:16.0 1st Qu.:0.0000
## Median :0.0000 Median :0.0000 Median :17.0 Median :0.0000
## Mean :0.1165 Mean :0.4734 Mean :16.7 Mean :0.2228
## 3rd Qu.:0.0000 3rd Qu.:1.0000 3rd Qu.:18.0 3rd Qu.:0.0000
## Max. :1.0000 Max. :1.0000 Max. :22.0 Max. :1.0000
## famsize Pstatus Mjob reason
## Min. :0.0000 Min. :0.0000 Min. :0.000 Min. :0.000
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:2.000 1st Qu.:0.000
## Median :1.0000 Median :0.0000 Median :2.499 Median :1.000
## Mean :0.7114 Mean :0.1038 Mean :2.499 Mean :1.273
## 3rd Qu.:1.0000 3rd Qu.:0.0000 3rd Qu.:4.000 3rd Qu.:2.000
## Max. :1.0000 Max. :1.0000 Max. :4.000 Max. :3.000
## traveltime studytime failures schoolsup
## Min. :1.000 Min. :1.000 Min. :0.0000 Min. :0.0000
## 1st Qu.:1.000 1st Qu.:2.000 1st Qu.:0.0000 1st Qu.:1.0000
## Median :1.000 Median :2.000 Median :0.0000 Median :1.0000
## Mean :1.448 Mean :2.023 Mean :0.3342 Mean :0.8709
## 3rd Qu.:2.000 3rd Qu.:2.023 3rd Qu.:0.0000 3rd Qu.:1.0000
## Max. :4.000 Max. :4.000 Max. :3.0000 Max. :1.0000
## famsup paid activities nursery
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000
## Median :0.0000 Median :1.0000 Median :0.0000 Median :0.0000
## Mean :0.3873 Mean :0.5418 Mean :0.4911 Mean :0.2051
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:0.0000
## Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000
## higher internet romantic famrel
## Min. :0.00000 Min. :0.0000 Min. :0.0000 Min. :1.000
## 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:4.000
## Median :0.00000 Median :0.0000 Median :1.0000 Median :4.000
## Mean :0.05063 Mean :0.1671 Mean :0.6658 Mean :3.944
## 3rd Qu.:0.00000 3rd Qu.:0.0000 3rd Qu.:1.0000 3rd Qu.:5.000
## Max. :1.00000 Max. :1.0000 Max. :1.0000 Max. :5.000
## freetime goout Dalc Walc
## Min. :1.000 Min. :1.000 Min. :1.000 Min. :1.000
## 1st Qu.:3.000 1st Qu.:2.000 1st Qu.:1.000 1st Qu.:1.000
## Median :3.000 Median :3.000 Median :1.000 Median :2.000
## Mean :3.235 Mean :3.109 Mean :1.481 Mean :2.291
## 3rd Qu.:4.000 3rd Qu.:4.000 3rd Qu.:2.000 3rd Qu.:3.000
## Max. :5.000 Max. :5.000 Max. :5.000 Max. :5.000
## health absences G3
## Min. :1.000 Min. : 0.000 Min. : 0.00
## 1st Qu.:3.000 1st Qu.: 0.000 1st Qu.: 8.00
## Median :4.000 Median : 4.000 Median :11.00
## Mean :3.554 Mean : 5.709 Mean :10.42
## 3rd Qu.:5.000 3rd Qu.: 8.000 3rd Qu.:14.00
## Max. :5.000 Max. :75.000 Max. :20.002.5 Outliers
Para detectar los valores atípicos (outliers), hay que hacer un análisis exploratorio gráfico previo, construyendo los boxplots de todas las variables.
# Boxplots of all variables together
# This visualization is not the best due to the difference between the scales
boxplot(data_xlsx, main="Outliers",
xlab="All explanatory variables",
ylab="Values",
col=c(1:ncol(data_xlsx)))
Se ha tomado la decisión de sustituir los outliers por la media de sus variables. Para ello, se ha construido la función outlier(), que realiza la detección y manipulación de los mismos. Cabe destacar que la decisión de sustituir los outliers por la media requeriría un análisis previo de la causa de estos valores atípicos.
# Recursive function that modifies outliers by the mean of their variable
outlier<-function(data, na.rm=T){
H<-1.5*IQR(data)
data[data<quantile(data, 0.25, na.rm = T)-H]<-NA
data[data>quantile(data, 0.75, na.rm = T)+H]<-NA
data[is.na(data)]<-mean(data, na.rm = T)
H<-1.5*IQR(data)
if (TRUE %in% (data<quantile(data, 0.25, na.rm = T)-H) |
TRUE %in% (data>quantile(data, 0.75, na.rm = T)+H))
outlier(data)
else
return(data)
}
# This data.frame is to preserve original data once the outliers are modified
original_data_xlsx<-data_xlsx
# Called to outlier function for each variable identified with outliers
data_xlsx<-as.data.frame(apply(data_xlsx, 2, outlier))
# We compare the original data and the fixed ones with respective boxplots
par(mfrow=c(1,2))
# Boxplot original data
boxplot(original_data_xlsx, main="Original data",
xlab="All explanatory variables",
ylab="Values",
col=c(1:ncol(original_data_xlsx)))
# Boxplot fixed data
boxplot(data_xlsx, main="Data with no outliers",
xlab="All explanatory variables",
ylab="Values",
col=c(1:ncol(data_xlsx)))
Vemos que se aprecian muy bien los boxplots debido a los altos valores de la penúltima variable (absences) en comparación con el resto. A continuación se repite el análisis de outliers pero con los datos normalizados para que se pueda tener una mejor visualización de las gráficas anteriores.
# Normalized original data
normalized_data_xlsx<-original_data_xlsx
normalized_data_xlsx<-scale(original_data_xlsx)
# Boxplots of all variables together
# This visualization is not the best due to the difference between the scales
boxplot(normalized_data_xlsx, main="Outliers",
xlab="All explanatory variables",
ylab="Values",
col=c(1:ncol(normalized_data_xlsx)))
# This data.frame is to preserve original data once the outliers are modified
normalized_original_data_xlsx<-normalized_data_xlsx
# Called to outlier function for each variable identified with outliers
normalized_data_xlsx<-as.data.frame(apply(normalized_data_xlsx, 2, outlier))
# We compare the normalized original data and the fixed ones with respective boxplots
par(mfrow=c(1,2))
# Boxplot normalized original data
boxplot(normalized_original_data_xlsx, main="Normalized original data",
xlab="All explanatory variables",
ylab="Values",
col=c(1:ncol(normalized_original_data_xlsx)))
# Boxplot fixed data
boxplot(normalized_data_xlsx, main="Normalized data with no outliers",
xlab="All explanatory variables",
ylab="Values",
col=c(1:ncol(normalized_data_xlsx)))
3 Análisis de Componentes Principales (ACP)
3.1 Correlación
En primer lugar es necesario comprobar que las variables no son independientes. A nivel de la muestra recogida en la base de datos, esto se puede hacer calculando y observando la matriz de correlaciones. A nivel poblacional se justificará que hay correlación realizando el test de Bartlett (el contraste de esfericidad de Bartlett permite comprobar si las correlaciones son distintas de 0 de modo significativo. La hipótesis nula es que \(det(R)=1\) , siendo \(R\) la matriz de correlaciones). El siguiente código realiza las dos comprobaciones indicadas.
###############################
# Correlation at sample level #
###############################
# Are the variables correlated at sample level?
correlation_matrix<-cor(data_xlsx)
correlation_matrix## school sex age address famsize Pstatus
## school 1 NA NA NA NA NA
## sex NA 1.000000000 -0.04064933 NA -0.0898618018 NA
## age NA -0.040649328 1.00000000 NA -0.0455885843 NA
## address NA NA NA 1 NA NA
## famsize NA -0.089861802 -0.04558858 NA 1.0000000000 NA
## Pstatus NA NA NA NA NA 1
## Mjob NA -0.102773909 0.08962071 NA 0.0915722328 NA
## reason NA 0.009821817 -0.01340060 NA 0.0009812261 NA
## traveltime NA 0.020908545 0.11516528 NA -0.0577081482 NA
## studytime NA 0.187170093 -0.05072535 NA 0.0341291336 NA
## failures NA NA NA NA NA NA
## schoolsup NA NA NA NA NA NA
## famsup NA 0.151622617 0.13021861 NA -0.1128930671 NA
## paid NA 0.129125619 0.02681464 NA -0.0138821617 NA
## activities NA -0.099833177 0.09440942 NA -0.0001131763 NA
## nursery NA NA NA NA NA NA
## higher NA NA NA NA NA NA
## internet NA NA NA NA NA NA
## romantic NA 0.102023009 -0.15316341 NA 0.0343948162 NA
## famrel NA 0.069375471 0.03237037 NA 0.0069223355 NA
## freetime NA 0.207351483 0.02663625 NA -0.0131216979 NA
## goout NA 0.075897396 0.11147582 NA -0.0230643843 NA
## Dalc NA 0.224996943 0.09032711 NA -0.1485945736 NA
## Walc NA 0.274193767 0.09719130 NA -0.1034250076 NA
## health NA 0.143588173 -0.04372758 NA 0.0289916591 NA
## absences NA 0.006282370 0.13691739 NA -0.1095387461 NA
## G3 NA 0.103455647 -0.15955029 NA -0.0814071098 NA
## Mjob reason traveltime studytime failures
## school NA NA NA NA NA
## sex -0.10277391 0.0098218168 0.020908545 0.18717009 NA
## age 0.08962071 -0.0134006027 0.115165279 -0.05072535 NA
## address NA NA NA NA NA
## famsize 0.09157223 0.0009812261 -0.057708148 0.03412913 NA
## Pstatus NA NA NA NA NA
## Mjob 1.00000000 -0.0486185847 0.105091345 0.03289251 NA
## reason -0.04861858 1.0000000000 0.083587795 0.10945933 NA
## traveltime 0.10509135 0.0835877947 1.000000000 0.01494850 NA
## studytime 0.03289251 0.1094593263 0.014948499 1.00000000 NA
## failures NA NA NA NA 1
## schoolsup NA NA NA NA NA
## famsup 0.14047224 0.0762296929 0.006122867 0.08144172 NA
## paid 0.16676148 0.0867430657 0.040906720 0.12737933 NA
## activities 0.12072516 0.0207472217 0.007557387 0.06199645 NA
## nursery NA NA NA NA NA
## higher NA NA NA NA NA
## internet NA NA NA NA NA
## romantic -0.04106315 -0.0050504430 0.029087153 0.13600710 NA
## famrel 0.02556143 0.0572004768 -0.024250888 -0.04504248 NA
## freetime -0.05205738 0.0352369814 -0.046113188 0.07137651 NA
## goout -0.01377800 0.0170752616 -0.037824830 0.01081643 NA
## Dalc -0.05968205 0.0297294798 0.038855019 0.11420202 NA
## Walc -0.01954561 0.0602322517 0.030763971 0.12417235 NA
## health -0.04147519 0.0700755406 -0.019762406 0.05172415 NA
## absences -0.04718168 -0.0302451714 -0.052442404 -0.02821618 NA
## G3 -0.14572680 -0.0085015237 -0.105385244 0.03982203 NA
## schoolsup famsup paid activities nursery higher
## school NA NA NA NA NA NA
## sex NA 0.151622617 0.129125619 -0.0998331770 NA NA
## age NA 0.130218615 0.026814641 0.0944094170 NA NA
## address NA NA NA NA NA NA
## famsize NA -0.112893067 -0.013882162 -0.0001131763 NA NA
## Pstatus NA NA NA NA NA NA
## Mjob NA 0.140472242 0.166761480 0.1207251610 NA NA
## reason NA 0.076229693 0.086743066 0.0207472217 NA NA
## traveltime NA 0.006122867 0.040906720 0.0075573873 NA NA
## studytime NA 0.081441723 0.127379327 0.0619964523 NA NA
## failures NA NA NA NA NA NA
## schoolsup 1 NA NA NA NA NA
## famsup NA 1.000000000 0.293184339 -0.0015001081 NA NA
## paid NA 0.293184339 1.000000000 -0.0213823757 NA NA
## activities NA -0.001500108 -0.021382376 1.0000000000 NA NA
## nursery NA NA NA NA 1 NA
## higher NA NA NA NA NA 1
## internet NA NA NA NA NA NA
## romantic NA 0.012439892 0.005535859 0.0196505434 NA NA
## famrel NA 0.025926895 0.036491388 -0.0350101150 NA NA
## freetime NA 0.025196125 0.070312152 -0.0532456586 NA NA
## goout NA 0.015631443 -0.010493268 -0.0460876858 NA NA
## Dalc NA 0.083896122 0.013793327 0.0777486766 NA NA
## Walc NA 0.086687935 -0.060453643 0.0374766962 NA NA
## health NA -0.029297452 0.078132403 -0.0239228148 NA NA
## absences NA -0.024280907 -0.058049073 -0.0526057549 NA NA
## G3 NA 0.039157145 -0.101996241 -0.0160997013 NA NA
## internet romantic famrel freetime goout
## school NA NA NA NA NA
## sex NA 0.102023009 0.069375471 0.207351483 0.075897396
## age NA -0.153163415 0.032370369 0.026636252 0.111475815
## address NA NA NA NA NA
## famsize NA 0.034394816 0.006922336 -0.013121698 -0.023064384
## Pstatus NA NA NA NA NA
## Mjob NA -0.041063155 0.025561429 -0.052057378 -0.013777997
## reason NA -0.005050443 0.057200477 0.035236981 0.017075262
## traveltime NA 0.029087153 -0.024250888 -0.046113188 -0.037824830
## studytime NA 0.136007100 -0.045042484 0.071376514 0.010816431
## failures NA NA NA NA NA
## schoolsup NA NA NA NA NA
## famsup NA 0.012439892 0.025926895 0.025196125 0.015631443
## paid NA 0.005535859 0.036491388 0.070312152 -0.010493268
## activities NA 0.019650543 -0.035010115 -0.053245659 -0.046087686
## nursery NA NA NA NA NA
## higher NA NA NA NA NA
## internet 1 NA NA NA NA
## romantic NA 1.000000000 0.045505206 0.032679335 -0.007869926
## famrel NA 0.045505206 1.000000000 0.093644457 0.087047617
## freetime NA 0.032679335 0.093644457 1.000000000 0.236439409
## goout NA -0.007869926 0.087047617 0.236439409 1.000000000
## Dalc NA 0.067951491 -0.065351965 0.130440135 0.210311226
## Walc NA 0.010140952 -0.093566930 0.125357672 0.420385745
## health NA -0.026342317 0.043777985 0.065140450 -0.009577254
## absences NA -0.059041770 -0.023375472 0.015057208 0.131320084
## G3 NA 0.129969950 0.057689643 -0.003327928 -0.132791474
## Dalc Walc health absences G3
## school NA NA NA NA NA
## sex 0.22499694 0.27419377 0.143588173 0.00628237 0.103455647
## age 0.09032711 0.09719130 -0.043727578 0.13691739 -0.159550295
## address NA NA NA NA NA
## famsize -0.14859457 -0.10342501 0.028991659 -0.10953875 -0.081407110
## Pstatus NA NA NA NA NA
## Mjob -0.05968205 -0.01954561 -0.041475191 -0.04718168 -0.145726800
## reason 0.02972948 0.06023225 0.070075541 -0.03024517 -0.008501524
## traveltime 0.03885502 0.03076397 -0.019762406 -0.05244240 -0.105385244
## studytime 0.11420202 0.12417235 0.051724153 -0.02821618 0.039822028
## failures NA NA NA NA NA
## schoolsup NA NA NA NA NA
## famsup 0.08389612 0.08668793 -0.029297452 -0.02428091 0.039157145
## paid 0.01379333 -0.06045364 0.078132403 -0.05804907 -0.101996241
## activities 0.07774868 0.03747670 -0.023922815 -0.05260575 -0.016099701
## nursery NA NA NA NA NA
## higher NA NA NA NA NA
## internet NA NA NA NA NA
## romantic 0.06795149 0.01014095 -0.026342317 -0.05904177 0.129969950
## famrel -0.06535196 -0.09356693 0.043777985 -0.02337547 0.057689643
## freetime 0.13044013 0.12535767 0.065140450 0.01505721 -0.003327928
## goout 0.21031123 0.42038575 -0.009577254 0.13132008 -0.132791474
## Dalc 1.00000000 0.53544088 0.074332685 0.05337669 -0.078730632
## Walc 0.53544088 1.00000000 0.092476317 0.17807684 -0.051939324
## health 0.07433269 0.09247632 1.000000000 -0.04993949 -0.061334605
## absences 0.05337669 0.17807684 -0.049939495 1.00000000 0.121485428
## G3 -0.07873063 -0.05193932 -0.061334605 0.12148543 1.000000000Vemos que hay algunas filas en la matriz de correlación con el valor ‘NA’ correspondientes a las variables school, address, Pstatus, failures, schoolsup, nursery, higher e internet. Esto se debe a que todos los registros de estas variables tienen el mismo valor (en este caso esto se debe al tratamiento de los outliers), haciendo que su desviación típica sea igual a 0 y, por lo tanto, no podemos calcular su correlación con otras variables (ya que sería dividir por 0, lo cual no es posible). Por tanto, debido a que disponemos de muchas más variables, se ha decidido eliminar estas variables.
# The variables 'school', 'address', 'Pstatus', 'failures', 'schoolsup', 'nursery', 'higher' and 'internet' are eliminated
data_xlsx_reduced<-data_xlsx[, -c(1, 4, 6, 11, 12, 16, 17, 18)]
# The first three records in the database are displayed
head(data_xlsx_reduced, n=10)En cualquier caso, para hacer el análisis de la correlación hemos trabajado con la matriz de correlaciones de todas las variables. Aquí hay un problema porque estamos mezclando variables por un lado continuas o discretas, con variables categóricas que han sido codificadas (sexo, school, address, etc.). Es un problema porque por defecto estamos obteniendo las correlaciones de Pearson que trata todas las variables como cuantitativas y no lo son, llevando así a engaño. Para variables categóricas habría que cambiar el parámetro de la función de R para que obtuviera la correlación de Spearman. Pero claro, ya no podríamos hacer una reducción de la dimensión con todas a la vez. Por un lado lo haríamos con las cuantitativas que es para lo que funciona el ACP por su construcción y por otro lado investigaríamos como hacerlo con variables categóricas (lo cual no hemos visto). Este comentario es solo para tenerlo en cuenta de cara al futuro. Para esta práctica, el procedimiento seguido mezclando todas las variables es suficiente. En conclusión, recordar para trabajos futuros que la reducción de la dimensión la hacéis con las variables del mismo tipo (cuantitativas).
De acuerdo con los resultados numéricos a continuación, se observa que los datos no están correlacionados ni a nivel muestral (ver matriz de correlación) ni a nivel poblacional (la prueba de esfericidad de Bartlett no es significativa).
###############################
# Correlation at sample level #
###############################
# Are the variables correlated at sample level?
correlation_matrix<-cor(data_xlsx_reduced)
correlation_matrix## sex age famsize Mjob reason
## sex 1.000000000 -0.04064933 -0.0898618018 -0.10277391 0.0098218168
## age -0.040649328 1.00000000 -0.0455885843 0.08962071 -0.0134006027
## famsize -0.089861802 -0.04558858 1.0000000000 0.09157223 0.0009812261
## Mjob -0.102773909 0.08962071 0.0915722328 1.00000000 -0.0486185847
## reason 0.009821817 -0.01340060 0.0009812261 -0.04861858 1.0000000000
## traveltime 0.020908545 0.11516528 -0.0577081482 0.10509135 0.0835877947
## studytime 0.187170093 -0.05072535 0.0341291336 0.03289251 0.1094593263
## famsup 0.151622617 0.13021861 -0.1128930671 0.14047224 0.0762296929
## paid 0.129125619 0.02681464 -0.0138821617 0.16676148 0.0867430657
## activities -0.099833177 0.09440942 -0.0001131763 0.12072516 0.0207472217
## romantic 0.102023009 -0.15316341 0.0343948162 -0.04106315 -0.0050504430
## famrel 0.069375471 0.03237037 0.0069223355 0.02556143 0.0572004768
## freetime 0.207351483 0.02663625 -0.0131216979 -0.05205738 0.0352369814
## goout 0.075897396 0.11147582 -0.0230643843 -0.01377800 0.0170752616
## Dalc 0.224996943 0.09032711 -0.1485945736 -0.05968205 0.0297294798
## Walc 0.274193767 0.09719130 -0.1034250076 -0.01954561 0.0602322517
## health 0.143588173 -0.04372758 0.0289916591 -0.04147519 0.0700755406
## absences 0.006282370 0.13691739 -0.1095387461 -0.04718168 -0.0302451714
## G3 0.103455647 -0.15955029 -0.0814071098 -0.14572680 -0.0085015237
## traveltime studytime famsup paid activities
## sex 0.020908545 0.18717009 0.151622617 0.129125619 -0.0998331770
## age 0.115165279 -0.05072535 0.130218615 0.026814641 0.0944094170
## famsize -0.057708148 0.03412913 -0.112893067 -0.013882162 -0.0001131763
## Mjob 0.105091345 0.03289251 0.140472242 0.166761480 0.1207251610
## reason 0.083587795 0.10945933 0.076229693 0.086743066 0.0207472217
## traveltime 1.000000000 0.01494850 0.006122867 0.040906720 0.0075573873
## studytime 0.014948499 1.00000000 0.081441723 0.127379327 0.0619964523
## famsup 0.006122867 0.08144172 1.000000000 0.293184339 -0.0015001081
## paid 0.040906720 0.12737933 0.293184339 1.000000000 -0.0213823757
## activities 0.007557387 0.06199645 -0.001500108 -0.021382376 1.0000000000
## romantic 0.029087153 0.13600710 0.012439892 0.005535859 0.0196505434
## famrel -0.024250888 -0.04504248 0.025926895 0.036491388 -0.0350101150
## freetime -0.046113188 0.07137651 0.025196125 0.070312152 -0.0532456586
## goout -0.037824830 0.01081643 0.015631443 -0.010493268 -0.0460876858
## Dalc 0.038855019 0.11420202 0.083896122 0.013793327 0.0777486766
## Walc 0.030763971 0.12417235 0.086687935 -0.060453643 0.0374766962
## health -0.019762406 0.05172415 -0.029297452 0.078132403 -0.0239228148
## absences -0.052442404 -0.02821618 -0.024280907 -0.058049073 -0.0526057549
## G3 -0.105385244 0.03982203 0.039157145 -0.101996241 -0.0160997013
## romantic famrel freetime goout Dalc
## sex 0.102023009 0.069375471 0.207351483 0.075897396 0.22499694
## age -0.153163415 0.032370369 0.026636252 0.111475815 0.09032711
## famsize 0.034394816 0.006922336 -0.013121698 -0.023064384 -0.14859457
## Mjob -0.041063155 0.025561429 -0.052057378 -0.013777997 -0.05968205
## reason -0.005050443 0.057200477 0.035236981 0.017075262 0.02972948
## traveltime 0.029087153 -0.024250888 -0.046113188 -0.037824830 0.03885502
## studytime 0.136007100 -0.045042484 0.071376514 0.010816431 0.11420202
## famsup 0.012439892 0.025926895 0.025196125 0.015631443 0.08389612
## paid 0.005535859 0.036491388 0.070312152 -0.010493268 0.01379333
## activities 0.019650543 -0.035010115 -0.053245659 -0.046087686 0.07774868
## romantic 1.000000000 0.045505206 0.032679335 -0.007869926 0.06795149
## famrel 0.045505206 1.000000000 0.093644457 0.087047617 -0.06535196
## freetime 0.032679335 0.093644457 1.000000000 0.236439409 0.13044013
## goout -0.007869926 0.087047617 0.236439409 1.000000000 0.21031123
## Dalc 0.067951491 -0.065351965 0.130440135 0.210311226 1.00000000
## Walc 0.010140952 -0.093566930 0.125357672 0.420385745 0.53544088
## health -0.026342317 0.043777985 0.065140450 -0.009577254 0.07433269
## absences -0.059041770 -0.023375472 0.015057208 0.131320084 0.05337669
## G3 0.129969950 0.057689643 -0.003327928 -0.132791474 -0.07873063
## Walc health absences G3
## sex 0.27419377 0.143588173 0.00628237 0.103455647
## age 0.09719130 -0.043727578 0.13691739 -0.159550295
## famsize -0.10342501 0.028991659 -0.10953875 -0.081407110
## Mjob -0.01954561 -0.041475191 -0.04718168 -0.145726800
## reason 0.06023225 0.070075541 -0.03024517 -0.008501524
## traveltime 0.03076397 -0.019762406 -0.05244240 -0.105385244
## studytime 0.12417235 0.051724153 -0.02821618 0.039822028
## famsup 0.08668793 -0.029297452 -0.02428091 0.039157145
## paid -0.06045364 0.078132403 -0.05804907 -0.101996241
## activities 0.03747670 -0.023922815 -0.05260575 -0.016099701
## romantic 0.01014095 -0.026342317 -0.05904177 0.129969950
## famrel -0.09356693 0.043777985 -0.02337547 0.057689643
## freetime 0.12535767 0.065140450 0.01505721 -0.003327928
## goout 0.42038575 -0.009577254 0.13132008 -0.132791474
## Dalc 0.53544088 0.074332685 0.05337669 -0.078730632
## Walc 1.00000000 0.092476317 0.17807684 -0.051939324
## health 0.09247632 1.000000000 -0.04993949 -0.061334605
## absences 0.17807684 -0.049939495 1.00000000 0.121485428
## G3 -0.05193932 -0.061334605 0.12148543 1.000000000det(correlation_matrix)## [1] 0.1963095###################################
# Correlation at population level #
###################################
# Bartlett's sphericity test:
# This test checks whether the correlations are significantly different from 0
# The null hypothesis is H_0; det(R)=1 means the variables are uncorrelated
# R denotes the correlation matrix
# cortest.bartlett function in the package pysch performs this test
# This function works with standardized data.
# Standardization
data_xlsx_scale<-scale(data_xlsx_reduced)
# Bartlett's sphericity test
cortest.bartlett(cor(data_xlsx_scale))## $chisq
## [1] 149.5104
##
## $p.value
## [1] 0.880634
##
## $df
## [1] 171A la vista de los resultados del test (p-value > 0.001), no se puede rechazar que los datos sean incorrelados, es decir, que sean independientes. Por tanto, no tendría mucho sentido realmente plantearse una reducción de la dimensión mediante el procedimiento de Análisis de Componentes Principales (ACP) o Análisis Factorial (AF). Sin embargo, vamos a realizarlo con el fin de mostrar cómo se haría.
3.2 Realización del ACP
El siguiente código realiza el ACP, obteniendo los vectores propios que generan cada componente, así como sus valores propios, que corresponden a la varianza de cada una.
# The 'prcomp' function in the base R package performs this analysis
# Parameters 'scale' and 'center' are set to TRUE to consider standardized data
PCA<-prcomp(data_xlsx_reduced, scale=T, center=T)
# The field 'rotation' of the 'PCA' object is a matrix
# Its columns are the coefficients of the principal components
# Indicates the weight of each variable in the corresponding principal component
PCA$rotation## PC1 PC2 PC3 PC4 PC5
## sex 0.357988987 -0.138411478 0.322697872 -0.054267662 0.067068174
## age 0.131267998 0.339957291 -0.331365990 -0.176090198 0.145042062
## famsize -0.161590831 0.071228092 0.069841422 0.498471643 0.008965111
## Mjob -0.044001488 0.492619024 0.007931548 -0.002584637 -0.034385460
## reason 0.094556511 0.101828074 0.186479728 0.076449751 -0.014264218
## traveltime 0.038760178 0.271918073 -0.015201624 -0.029883779 -0.201107849
## studytime 0.194298544 0.043651596 0.346186243 0.055168073 -0.299681437
## famsup 0.182311516 0.318048173 0.250930379 -0.422734663 0.118971288
## paid 0.104874161 0.392601699 0.366688720 -0.111221050 0.194371680
## activities -0.001007259 0.198060733 -0.082199948 -0.009605777 -0.425364174
## romantic 0.054747386 -0.182091711 0.321552812 0.017569240 -0.268945880
## famrel 0.008972162 -0.007201713 0.152799749 0.037058002 0.483625472
## freetime 0.270209598 -0.087769965 0.102905512 0.197956338 0.343553056
## goout 0.366548664 -0.025734817 -0.255907511 0.205602698 0.234858149
## Dalc 0.465023159 -0.022189337 -0.104584992 0.036943252 -0.268693347
## Walc 0.519812584 -0.059801406 -0.197971823 0.062367826 -0.178380077
## health 0.123665610 -0.012696051 0.178673612 0.349364636 0.083594720
## absences 0.139648827 -0.166723530 -0.298530527 -0.351770959 0.135539274
## G3 -0.049366010 -0.403967956 0.219294237 -0.427111617 -0.032997838
## PC6 PC7 PC8 PC9 PC10
## sex -0.10391169 -0.10839454 0.16365183 -0.139257379 0.20557884
## age 0.02303033 0.15411088 -0.07849430 -0.088680141 0.22237639
## famsize 0.31314158 -0.10650027 -0.10643092 0.301916325 0.25727557
## Mjob 0.33116460 -0.11367398 0.06714274 -0.011128482 0.24163018
## reason -0.28271349 0.52853728 -0.45083302 0.395241591 -0.25834592
## traveltime -0.28597132 0.48146496 0.49939418 0.102358812 0.29838847
## studytime 0.14566345 -0.03532613 -0.19240446 0.298003384 0.20817459
## famsup 0.04890066 -0.15759837 -0.04364686 0.004294231 -0.21485205
## paid -0.04884339 -0.22592945 -0.01539739 0.099576591 -0.05774483
## activities 0.29229668 0.16916938 -0.43267571 -0.456748314 -0.06963954
## romantic 0.37116072 0.22067039 0.32392297 -0.002982180 0.04013603
## famrel 0.23783650 0.48463424 -0.04143655 -0.345570359 0.13434361
## freetime 0.18335791 0.05011800 0.04997650 -0.026074757 -0.15056639
## goout 0.26727354 0.06983414 0.05217576 0.171907752 -0.17261233
## Dalc -0.05127171 -0.03837628 0.05303724 -0.175134786 -0.11882933
## Walc 0.01527602 -0.04995084 -0.02075521 0.044733675 0.02095125
## health -0.41945363 -0.13236526 -0.27465340 -0.343906254 0.42252253
## absences 0.07709932 -0.02275517 -0.22455581 0.325558847 0.48661653
## G3 0.15387047 0.10182366 -0.17674737 -0.028642312 0.16806270
## PC11 PC12 PC13 PC14 PC15
## sex -0.143978916 0.197058675 -0.263363032 -0.002061750 -0.30748633
## age -0.324515844 0.487108748 0.246267697 -0.055017486 0.22599694
## famsize 0.035308463 0.488004212 -0.100712430 0.301302539 -0.24592511
## Mjob 0.242299180 -0.380084492 -0.340510269 0.140047267 0.23154673
## reason 0.136849347 -0.006144228 -0.028546872 0.161250872 0.01531048
## traveltime -0.128723700 -0.096826268 -0.131715638 0.175449808 -0.08922372
## studytime -0.306801210 -0.053734447 -0.089034545 -0.600946157 0.23786877
## famsup 0.169697930 0.312720452 -0.037366643 0.189046850 0.28369633
## paid 0.013212497 -0.190954934 0.330258017 -0.002415829 -0.48699543
## activities -0.195063366 -0.123660716 0.051135496 0.131056496 -0.27602362
## romantic 0.205147170 0.046890289 0.599315081 0.119626838 0.17559823
## famrel 0.253285005 0.053805328 -0.121153360 -0.360065071 -0.15430665
## freetime -0.608901579 -0.273425928 0.012434492 0.362500709 0.14683036
## goout 0.173901756 -0.130000373 0.008027874 -0.067047736 0.10288304
## Dalc 0.133241456 0.105844736 0.006299095 0.002167397 -0.16751420
## Walc 0.240014167 0.043755572 -0.175019668 0.064239619 -0.01442977
## health 0.177947027 -0.109242105 0.214537433 0.158185734 0.33059549
## absences 0.057199275 -0.238459491 0.267427336 0.061761636 -0.21955060
## G3 -0.004378558 0.027574136 -0.287120761 0.320784978 0.09547792
## PC16 PC17 PC18 PC19
## sex 0.04393326 -0.59012512 0.15568682 0.17953965
## age 0.08134033 0.01016004 0.39072602 -0.02951198
## famsize -0.03244157 0.13339112 -0.12771368 0.03680353
## Mjob 0.32382919 -0.11862328 0.22129067 0.09097057
## reason 0.21928319 -0.18524533 0.14898620 0.06297112
## traveltime -0.27960227 0.15174886 -0.19130319 0.02558064
## studytime -0.07042271 0.13531241 -0.09935796 0.01473493
## famsup -0.15542534 -0.04212068 -0.51168242 0.05009488
## paid -0.15142989 0.26956896 0.29247492 -0.13899018
## activities -0.24644475 -0.17754971 -0.12756713 0.02002253
## romantic 0.10271605 -0.16220615 0.08151722 -0.05188973
## famrel 0.10093130 0.14617421 -0.17382129 -0.11696431
## freetime 0.20357093 0.09979195 -0.16600066 -0.10908413
## goout -0.57289816 -0.07682474 0.17385188 0.37020706
## Dalc 0.40773927 0.46952353 -0.04379517 0.44874450
## Walc -0.03853755 0.01675788 0.04086588 -0.74013052
## health -0.16081574 0.05391938 -0.03258438 0.06412728
## absences 0.13395249 -0.11290448 -0.31340469 0.09598673
## G3 -0.21108835 0.37180623 0.35777495 0.04881103# Standard deviations of each principal component
PCA$sdev## [1] 1.4966008 1.2622886 1.2484457 1.1234160 1.1212370 1.0524091 1.0247914
## [8] 1.0062424 0.9838951 0.9547318 0.9281544 0.8978759 0.8823035 0.8702234
## [15] 0.8226476 0.8007155 0.7737687 0.7478517 0.5998485Cada componente principal se obtiene de forma sencilla como una combinación lineal de todas las variables con los coeficientes indicados por las columnas de la matriz de rotación.
3.3 Tasa de variación explicada
# The function 'summary' applied to the 'PCA' object provides relevant information
# - Standard deviations of each principal component
# - Proportion of variance explained and cummulative variance
summary(PCA)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.4966 1.26229 1.24845 1.12342 1.12124 1.05241 1.02479
## Proportion of Variance 0.1179 0.08386 0.08203 0.06642 0.06617 0.05829 0.05527
## Cumulative Proportion 0.1179 0.20175 0.28378 0.35020 0.41637 0.47466 0.52994
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 1.00624 0.98390 0.95473 0.92815 0.89788 0.88230 0.87022
## Proportion of Variance 0.05329 0.05095 0.04797 0.04534 0.04243 0.04097 0.03986
## Cumulative Proportion 0.58323 0.63418 0.68215 0.72749 0.76992 0.81089 0.85075
## PC15 PC16 PC17 PC18 PC19
## Standard deviation 0.82265 0.80072 0.77377 0.74785 0.59985
## Proportion of Variance 0.03562 0.03374 0.03151 0.02944 0.01894
## Cumulative Proportion 0.88637 0.92011 0.95163 0.98106 1.00000Los siguientes gráficos ilustran el comportamiento de la varianza explicada por cada componente principal, así como el comportamiento de la varianza explicada acumulada.
# The following graph shows the proportion of explained variance
Explained_variance <- PCA$sdev^2 / sum(PCA$sdev^2)
p13<-ggplot(data = data.frame(Explained_variance, pc = 1:ncol(data_xlsx_reduced)),
aes(x = pc, y = Explained_variance, fill=Explained_variance )) +
geom_col(width = 0.3) + scale_y_continuous(limits = c(0,0.6)) + theme_bw() +
labs(x = "Principal component", y= "Proportion of variance")
# The following graph shows the proportion of cumulative explained variance
Cummulative_variance<-cumsum(Explained_variance)
p14<-ggplot( data = data.frame(Cummulative_variance, pc = 1:ncol(data_xlsx_reduced)),
aes(x = pc, y = Cummulative_variance ,fill=Cummulative_variance )) +
geom_col(width = 0.5) + scale_y_continuous(limits = c(0,1)) + theme_bw() +
labs(x = "Principal component",
y = "Proportion of cumulative variance")
p13
p14
3.4 Número de componentes principales apropiado
Para elegir el número adecuado de componentes principales, existen diferentes métodos. Nosotros vamos a utilizar la Regla de Abdi et al. (2010): Se promedian las varianzas explicadas por las componentes principales y se seleccionan aquellas cuya proporción de varianza explicada supera la media.
#######################
# Rule of Abdi et al. #
#######################
# Variances
PCA$sdev^2## [1] 2.2398140 1.5933726 1.5586166 1.2620636 1.2571724 1.1075649 1.0501974
## [8] 1.0125238 0.9680496 0.9115128 0.8614706 0.8061811 0.7784594 0.7572887
## [15] 0.6767490 0.6411453 0.5987180 0.5592821 0.3598182# Average of variances
mean(PCA$sdev^2)## [1] 1A la vista de los resultados, en este caso, se eligen ocho direcciones principales que, tal y como se puede ver, acumulan alrededor del 60% de varianza explicada.
3.5 Salidas gráficas de interés
Los siguientes gráficos muestran la comparativa entre distintas componentes principales mediante una proyección al plano sobre cada dos componentes. En esta representación se aprecia cuáles de las variables originales tienen mayor o menor peso en cada una de las componentes enfrentadas.
# These graphical outputs show the projection of the variables in two dimensions
# Display the weight of the variable in the direction of the principal component
p15<-fviz_pca_var(PCA, repel=TRUE, col.var="cos2",
legend.title="Distance", title="Variables")+theme_bw()
p16<-fviz_pca_var(PCA, axes=c(1,3), repel=TRUE, col.var="cos2",
legend.title="Distance", title="Variables")+theme_bw()
p17<-fviz_pca_var(PCA, axes=c(2,3), repel=TRUE, col.var="cos2",
legend.title="Distance", title="Variables")+theme_bw()
# Displaying graphics
p15
p16
p17
Es posible también representar las observaciones de los objetos junto con las componentes principales mediante la orden contrib() de la función fviz_pca_ind() anterior, así como identificar con colores aquellas observaciones que mayor varianza explican.
# It is also possible to represent the observations
# As well as identify with colors those observations that explain the greatest
# variance of the principal components
p18<-fviz_pca_ind(PCA,col.ind = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel=TRUE,legend.title="Contrib.var", title="Records")+theme_bw()
p19<-fviz_pca_ind(PCA,axes=c(1,3),col.ind = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel=TRUE,legend.title="Contrib.var", title="Records")+theme_bw()
p20<-fviz_pca_ind(PCA,axes=c(2,3),col.ind = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel=TRUE,legend.title="Contrib.var", title="Records")+theme_bw()
# Displaying graphics
p18
p19
p20
Finalmente, se puede obtener una representación conjunta de variables y observaciones que relaciona visualmente las posibles relaciones entre las observaciones, las contribuciones de los individuos a las varianzas y el peso de las variables en cada componente principal.
# Joint representation of variables and observations
# Relates the possible relationships between the contributions of the records
# to the variances of the components and the weight of the variables in each
# principal component
p21<-fviz_pca(PCA,alpha.ind ="contrib", col.var = "cos2",
col.ind="seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"),
repel=TRUE, legend.title="Distancia")+theme_bw()
p22<-fviz_pca(PCA,axes=c(1,3),alpha.ind ="contrib",
col.var = "cos2",col.ind="seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"),
repel=TRUE, legend.title="Distancia")+theme_bw()
p23<-fviz_pca(PCA,axes=c(2,3),alpha.ind ="contrib",
col.var = "cos2",col.ind="seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"),
repel=TRUE, legend.title="Distancia")+theme_bw()
# Displaying graphics
p21
p22
p23
Finalmente, dado que el objeto de este estudio es reducir la dimensión del conjunto de datos, es posible obtener las coordenadas de los datos originales en el nuevo sistema de referencia. De hecho, se almacenan desde que usamos la función prcomp() para crear la variable PCA.
head(PCA$x)## PC1 PC2 PC3 PC4 PC5 PC6
## [1,] -0.5428656 2.2596133 -0.42483568 -0.3808055 0.1973962 0.5663425
## [2,] -1.4738393 0.9734234 -0.03380767 0.7280407 0.6325970 1.0126359
## [3,] 0.3584064 -0.4617297 -0.43632964 -1.6441159 -1.2831378 -0.2222366
## [4,] -2.7481266 -1.8532500 -0.57242625 0.3913111 -0.6674298 -1.6632617
## [5,] -1.8208294 -0.3271368 -0.59114102 0.8691652 -0.8624869 0.8546701
## [6,] -0.0997650 -2.5648834 0.38175788 -0.9531831 1.4412018 -0.5339585
## PC7 PC8 PC9 PC10 PC11 PC12
## [1,] 0.9718671 -0.04616334 0.597973268 -0.2713659 0.2946118 0.3563613
## [2,] 1.0458847 -0.73175801 -0.358337362 -0.4003332 0.6240971 -0.1919391
## [3,] 0.8209019 -1.31891932 -0.048872389 -1.2455403 1.6704745 -1.0877667
## [4,] -2.0681437 -0.21915098 -0.006630797 0.2687390 0.3567702 0.3884771
## [5,] -0.7398849 -0.20685083 -1.425692474 0.9659082 0.6450230 -0.5232906
## [6,] 0.5215257 0.11646362 -1.123298254 1.2918069 0.2256947 -0.9883929
## PC13 PC14 PC15 PC16 PC17 PC18
## [1,] 1.4946749 0.89213985 -0.2480674 -0.8491035 -0.4345776 -0.45757939
## [2,] 1.2684351 -0.26645691 -1.0759564 0.4760184 -0.2633952 0.33962636
## [3,] 0.4009926 0.80278031 0.2609793 1.6434488 -0.5685151 -1.55488057
## [4,] -0.6810208 0.17717504 0.4006045 -0.7550371 0.8415472 -0.06199315
## [5,] 0.3234729 0.71503870 0.5434042 0.4081917 -0.2519953 -0.30382251
## [6,] 0.2993773 0.06754642 0.4105511 1.1864022 -0.8217480 -0.16359497
## PC19
## [1,] 0.65503362
## [2,] -0.02653478
## [3,] 0.03444780
## [4,] 0.40193837
## [5,] -0.38369841
## [6,] -0.31557997Recordemos que no hemos podido rechazar que los datos sean incorrelados, es decir, que sean independientes. Por tanto, la reducción de la dimensión mediante el procedimiento de Análisis de Componentes Principales (ACP) o Análisis Factorial (AF) no tendría mucho sentido realmente.
4 Análisis Factorial (AF)
4.1 Correlación
Es necesario comprobar que las variables no son independientes. Ya vimos en el apartado 3, al comprobar la correlación de las variables para el Análisis de Componentes Principales, que las variables sí son independientes tras calcular la matriz de correlación y realizar el test de Bartlett. Esto nos dice, volvemos a repetirlo, qque la reducción de la dimensión mediante el procedimiento de Análisis de Componentes Principales (ACP) o Análisis Factorial (AF) no tendría mucho sentido realmente. Sin embargo, igual que hicimos en el Análisis de Componentes Principales, vamos a realizarlo con el fin de mostrar cómo se haría.
Caben destacar también los siguientes resultados gráficos que proporcionan una idea intuitiva de la correlación entre las variables:
# Polychoric correlation matrix
# Remember that package 'polycor' is required for 'hetcor' function
poly_cor<-hetcor(data_xlsx_reduced)$correlations
# Remember that package 'ggcoorplot' is required for 'ggcorrplot' function
ggcorrplot(poly_cor, type="lower", hc.order=T)
# Another interesting visual representation is the following
# Remember that package 'corrplot' is required for 'corrplot' function
corrplot(cor(data_xlsx_reduced), order = "hclust", tl.col='black', tl.cex=1)
# Or this other option is also very visual
# Remember that package 'corrr' is required for 'rplot' function
data_xlsx_reduced_correlations <- correlate(data_xlsx_reduced)
rplot(data_xlsx_reduced_correlations, legend = TRUE, colours = c("firebrick1", "black","darkcyan"), print_cor = TRUE) 
Tal y como era de esperar, podemos observar que las variables apenas están correladas entre sí (lo cual corrobora los resultados numéricos de la matriz de correlación y el test de Bartlett que obtuvimos anteriormente).
4.2 Diferentes métodos de estimación
Para extraer los factores hay que escoger un método de estimación. Para ello, se utiliza la función fa(), que permite implementar hasta 6 métodos distintos.
A continuación, se comparan las salidas com el método del factor principal y con el método de máxima verosimilitud. Se calcula 1 factor, ya que es lo que se ve intuitivamente en los mapas de correlaciones anteriores.
### Test of two models with three factors
modelo1<-fa(poly_cor,
nfactors = 1,
rotate = "none",
fm="mle") # likelihood method
modelo2<-fa(poly_cor,
nfactors = 1,
rotate = "none",
fm="minres") # minimal residual model
# Outputs of these models: factorial matrices, etc.
print("Modelo 1: mle")## [1] "Modelo 1: mle"modelo1$loadings##
## Loadings:
## ML1
## sex 0.320
## age 0.123
## famsize -0.136
## Mjob
## reason
## traveltime
## studytime 0.151
## famsup 0.112
## paid
## activities
## romantic
## famrel
## freetime 0.185
## goout 0.455
## Dalc 0.603
## Walc 0.884
## health 0.105
## absences 0.185
## G3
##
## ML1
## SS loadings 1.623
## Proportion Var 0.085print("Modelo 2: minres")## [1] "Modelo 2: minres"modelo2$loadings##
## Loadings:
## MR1
## sex 0.376
## age 0.134
## famsize -0.163
## Mjob
## reason
## traveltime
## studytime 0.185
## famsup 0.161
## paid
## activities
## romantic
## famrel
## freetime 0.264
## goout 0.432
## Dalc 0.615
## Walc 0.809
## health 0.119
## absences 0.156
## G3
##
## MR1
## SS loadings 1.596
## Proportion Var 0.084Se comparan las comunalidades de ambos modelos para ver qué proporción de la varianza es explicada por los factores latentes.
# Comparing communalities
sort(modelo1$communality, decreasing = T)->c1
sort(modelo2$communality, decreasing = T)->c2
head(cbind(c1, c2))## c1 c2
## Walc 0.78204940 0.65515285
## Dalc 0.36417822 0.37848395
## goout 0.20704015 0.18671636
## sex 0.10217693 0.14110847
## freetime 0.03418673 0.06982291
## absences 0.03415778 0.03429621También se comparan las unicidades, esto es, la proporción de varianza que no ha sido explicada por el factor (1-comunalidad):
# Comparison of uniqueness, that is, the proportion of variance that has not
# been explained by the factor (1-communality)
sort(modelo1$uniquenesses, decreasing = T)->u1
sort(modelo2$uniquenesses, decreasing = T)->u2
head(cbind(u1,u2), n=10)## u1 u2
## romantic 0.9993603 0.9998923
## paid 0.9993580 0.9997853
## traveltime 0.9988124 0.9986512
## activities 0.9987844 0.9976418
## Mjob 0.9985758 0.9975256
## G3 0.9958093 0.9967005
## reason 0.9954460 0.9956642
## famrel 0.9943776 0.9930515
## health 0.9890089 0.9859061
## famsup 0.9875302 0.9819405Se observa que con el método de máxima verosimilitud se obtienen valores un poco más altos, luego los factores latentes obtenidos con este método explican mejor la varianza de las variables observadas.
4.3 Número de factores latentes apropiado
Existen diferentes criterios, entre los que destacan el scree plot (Cattel 1966) y el análisis paralelo (Horn 1965).
El método del codo consiste en representar el gráfico de sedimentación, que se obtiene al representar en las ordenadas los valores propios en orden decreciente y en las abscias los números correspondiente a las componentes principales o factores latentes (según se esté haciendo un ACP o un AF). Uniendo los puntos se obtiene una figura que empieza con fuerte pendiente para luego seguir con una ligera inclinación. El método indica que hay que quedarse con el número de factores que haya hasta el “codo” de la gráfica, donde se pasa de la fuerte pendiente a la ligera inclinación.
Se llama gráfico de sedimentación (scree plot en inglés), ya que parece el perfil de una montaña con una gran pendiente hasta llegar a la base, que es una meseta donde se acumulan los guijarros caídos desde la cumbre (donde sedimentan).
Según los siguientes resultados gráficos, 1 se considera el número óptimo de factores latentes (análisis paralelo), coincidiendo con el primer gráfico scree plot.
# Scree plot
scree(poly_cor)
# Parallel analysis
fa.parallel(poly_cor, n.obs=100, fa="fa", fm="ml")
## Parallel analysis suggests that the number of factors = 1 and the number of components = NAOtra forma de estudiar el número óptimo de factores es con el siguiente test de hipótesis, que contrasta si el número de factores es suficiente o no:
# Remember that package 'stats' is required for 'factanal' function
# This function only performs the mle method
factanal(data_xlsx_reduced, factors=1, rotation="varimax")##
## Call:
## factanal(x = data_xlsx_reduced, factors = 1, rotation = "varimax")
##
## Uniquenesses:
## sex age famsize Mjob reason traveltime studytime
## 0.898 0.985 0.981 0.999 0.995 0.999 0.977
## famsup paid activities romantic famrel freetime goout
## 0.988 0.999 0.999 0.999 0.994 0.966 0.793
## Dalc Walc health absences G3
## 0.636 0.218 0.989 0.966 0.996
##
## Loadings:
## Factor1
## sex 0.320
## age 0.123
## famsize -0.136
## Mjob
## reason
## traveltime
## studytime 0.151
## famsup 0.112
## paid
## activities
## romantic
## famrel
## freetime 0.185
## goout 0.455
## Dalc 0.603
## Walc 0.884
## health 0.105
## absences 0.185
## G3
##
## Factor1
## SS loadings 1.623
## Proportion Var 0.085
##
## Test of the hypothesis that 1 factor is sufficient.
## The chi square statistic is 341.42 on 152 degrees of freedom.
## The p-value is 1.36e-16Así, se acepta la hipótesis de que 1 factor es suficiente.
4.4 Realización del AF
Finalmente, se estima el modelo factorial con 1 factor implementando una rotación de tipo varimax para buscar una interpretación más simple de la realidad.
modelo_varimax<-fa(poly_cor, nfactors=1, rotate="varimax", fa="mle")
# The rotated factorial matrix is shown
print(modelo_varimax$loadings, cut=0) ##
## Loadings:
## MR1
## sex 0.376
## age 0.134
## famsize -0.163
## Mjob -0.050
## reason 0.083
## traveltime 0.037
## studytime 0.185
## famsup 0.161
## paid 0.066
## activities 0.010
## romantic 0.049
## famrel -0.015
## freetime 0.264
## goout 0.432
## Dalc 0.615
## Walc 0.809
## health 0.119
## absences 0.156
## G3 -0.057
##
## MR1
## SS loadings 1.596
## Proportion Var 0.084Visualmente se podría hacer el esfuerzo de ver con qué variables correlaciona cada uno de los factores en la matriz factorial, pero esto es muy tedioso, por lo que se utiliza la siguiente representación:
fa.diagram(modelo_varimax)
En el diagrama se observa que el único factor latente que tenemos correlaciona con las variables Walc, Dalc, goout y sex, siendo un resultado muy lógico que era de esperar.
Vemos que solo tenemos un factor latente y que solo correlaciona cuatro variables. Esto se debe a que (lo volvemos a recordar por última vez) no hemos podido rechazar que los datos sean incorrelados, es decir, que sean independientes. Por tanto, la reducción de la dimensión mediante el procedimiento de Análisis de Componentes Principales (ACP) o Análisis Factorial (AF) no tendría mucho sentido realmente.
5 Análisis Discriminante Lineal (ADL)
La base de datos contiene información sobre la nota final del estudiante en la asignatura Matemáticas. Así, sería interesante dar un método de clasificación para la calificación según los demás indicadores.
Para ello, haremos que la variable G3 sea categórica, con dos categorías:
Aprobado si la nota dada en la variable G3 es mayor o igual que 10.
Suspenso si la nota dada en la variable G3 es menor que 10.
calification = c()
for(k in 1:nrow(data_xlsx)) {
if(data_xlsx_reduced$G3[k] >= 10){
calification[k] = "aprobado"
}
else
calification[k] = "suspenso"
}
calification<-as.factor(calification)
# 'G3' now is a cathegorical variable
data_xlsx_G3_cathegorical = data_xlsx_reduced
data_xlsx_G3_cathegorical$G3 = calification
data_xlsx_G3_cathegorical$G3## [1] suspenso suspenso aprobado aprobado aprobado aprobado aprobado suspenso
## [9] aprobado aprobado suspenso aprobado aprobado aprobado aprobado aprobado
## [17] aprobado aprobado suspenso aprobado aprobado aprobado aprobado aprobado
## [25] suspenso suspenso aprobado aprobado aprobado aprobado aprobado aprobado
## [33] aprobado aprobado aprobado suspenso aprobado aprobado aprobado aprobado
## [41] aprobado aprobado aprobado aprobado suspenso suspenso aprobado aprobado
## [49] aprobado suspenso aprobado aprobado aprobado aprobado aprobado aprobado
## [57] aprobado aprobado suspenso aprobado aprobado aprobado suspenso suspenso
## [65] aprobado aprobado aprobado suspenso suspenso aprobado aprobado aprobado
## [73] suspenso aprobado aprobado aprobado aprobado aprobado aprobado suspenso
## [81] aprobado aprobado suspenso aprobado aprobado suspenso suspenso aprobado
## [89] aprobado suspenso suspenso aprobado suspenso aprobado aprobado aprobado
## [97] aprobado aprobado aprobado suspenso suspenso aprobado aprobado suspenso
## [105] aprobado aprobado suspenso aprobado aprobado aprobado aprobado aprobado
## [113] aprobado aprobado suspenso aprobado aprobado aprobado suspenso aprobado
## [121] aprobado aprobado aprobado aprobado suspenso aprobado aprobado suspenso
## [129] suspenso aprobado suspenso suspenso aprobado aprobado suspenso suspenso
## [137] suspenso suspenso aprobado aprobado suspenso suspenso aprobado aprobado
## [145] suspenso aprobado suspenso aprobado suspenso aprobado suspenso aprobado
## [153] aprobado suspenso aprobado suspenso aprobado aprobado aprobado aprobado
## [161] suspenso suspenso suspenso aprobado suspenso aprobado aprobado aprobado
## [169] suspenso aprobado suspenso aprobado aprobado suspenso suspenso suspenso
## [177] aprobado suspenso suspenso aprobado suspenso aprobado aprobado suspenso
## [185] aprobado aprobado aprobado aprobado suspenso aprobado aprobado suspenso
## [193] suspenso aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [201] aprobado aprobado aprobado suspenso aprobado suspenso suspenso aprobado
## [209] aprobado suspenso suspenso aprobado aprobado suspenso aprobado aprobado
## [217] suspenso suspenso suspenso aprobado suspenso suspenso aprobado aprobado
## [225] aprobado suspenso aprobado aprobado suspenso aprobado aprobado aprobado
## [233] suspenso aprobado suspenso aprobado aprobado aprobado aprobado suspenso
## [241] aprobado aprobado suspenso aprobado suspenso aprobado aprobado suspenso
## [249] suspenso aprobado suspenso aprobado suspenso suspenso aprobado suspenso
## [257] aprobado aprobado aprobado suspenso aprobado suspenso aprobado suspenso
## [265] suspenso aprobado aprobado aprobado aprobado suspenso suspenso aprobado
## [273] aprobado aprobado aprobado aprobado suspenso suspenso suspenso aprobado
## [281] suspenso aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [289] aprobado aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [297] suspenso suspenso aprobado aprobado aprobado aprobado aprobado aprobado
## [305] aprobado aprobado aprobado suspenso aprobado aprobado suspenso aprobado
## [313] aprobado aprobado aprobado aprobado suspenso suspenso aprobado aprobado
## [321] aprobado suspenso aprobado aprobado aprobado aprobado aprobado aprobado
## [329] suspenso aprobado suspenso aprobado suspenso suspenso suspenso aprobado
## [337] aprobado suspenso aprobado aprobado aprobado suspenso aprobado suspenso
## [345] aprobado aprobado aprobado suspenso aprobado aprobado suspenso aprobado
## [353] suspenso suspenso aprobado suspenso aprobado aprobado aprobado aprobado
## [361] aprobado aprobado aprobado aprobado aprobado aprobado aprobado suspenso
## [369] aprobado aprobado suspenso aprobado aprobado suspenso aprobado aprobado
## [377] aprobado aprobado aprobado aprobado aprobado suspenso aprobado suspenso
## [385] suspenso aprobado suspenso suspenso suspenso suspenso suspenso aprobado
## [393] suspenso aprobado suspenso
## Levels: aprobado suspensoPor otro lado, sabemos que el AD se basa en el concepto de distancia de Mahalanobis y se utiliza comúnmente para clasificar casos en grupos. Para realizar un AD, necesitamos al menos una variable categórica que actúe como tu variable dependiente (esta es la variable data_xlsx_G3_cathegorical$G3 en nuestro caso) y varias variables independientes numéricas que se usarán para predecir el grupo de clasificación.
Las variables categóricas binarias también pueden ser incluidas como variables independientes en el AD si se codifican adecuadamente, aunque su interpretación puede ser menos clara que con las variables continuas. De esta forma, para las variables independientes, vamos a considerar las variables cuantitativas y cualitativas que están codificadas numéricamente.
Vamos a trabajar con una partición del dataset, como conjunto de entrenamiento (con aproximadamente el 80% de los registros), sobre el que realizaremos el análisis discriminante y otra partición, como conjunto test (con el 20% de los registros restantes), con el que realizaremos la validación del modelo. Usaremos el dataset data_xlsx_G3_cathegorical definido a partir de data_xlsx_reduced (en los histogramas de la exploración gráfica de datos que se realiza a continuación se justifica el uso de este en vez del original).
# Partitioning the dataset: training (80%) + test (20%)
train_index<-createDataPartition(data_xlsx_G3_cathegorical$G3, p=0.80)$Resample1
train_data<-data_xlsx_G3_cathegorical[train_index,]
test_data<-data_xlsx_G3_cathegorical[-train_index,]5.1 Exploración gráfica de los datos
En primer lugar exploramos como de bien (o mal) clasifica en las calificaciones cada una de las variables explicativas consideradas de forma independiente. Para ello, dibujamos los histogramas superpuestos. Si los histogramas se separan, la variable considerada sería un buen clasificador individual para la calificación.
# Variables list
variables <- c("sex", "age", "famsize", "Mjob", "reason", "traveltime",
"studytime", "famsup", "paid", "activities", "romantic",
"famrel", "freetime", "goout", "Dalc", "Walc", "health",
"absences")
# Graphics list
graphics <- list()
for (var in variables) {
p <- ggplot(data = train_data, aes_string(x = var, fill = "G3")) +
geom_bar(position = "identity", alpha = 0.5)
graphics[[var]] <- p
}
ggarrange(plotlist = graphics, nrow = 5, ncol = 4, common.legend = TRUE)
En este sentido, a primera vista no se ve ninguna variable que diferencia muy bien la calificación.
Si hubiéramos usado el dataset data_xlsx, hubiésemos visto que hay histogramas de algunas variables con una sola barra. Esto es porque dichas variables tienen el mismo valor para todos los registros (consecuencia del tratamiento de outliers), haciendo que su varianza sea cero (lo cual nos daría problemas para el Test de normalidad univariante de Shapiro-Wilks, el cual realizaremos posteriormente). En concreto, estas variables son school, address, Pstatus, failures, schoolsup, nursery, higher e internet. Por ello, como tenemos un gran número de variables, se ha decidido eliminarlas, tal y como hicimos en ACP y AF, resultando así el dataset que ya hemos definido anteriormente data_xlsx_reduced (el cual estamos utilizando).
A continuación, se explora qué pares de variables separan mejor entre calificaciones.
pairs(x = train_data[, c("sex", "age", "famsize", "Mjob", "reason", "traveltime",
"studytime", "famsup", "paid", "activities", "romantic",
"famrel", "freetime", "goout", "Dalc", "Walc", "health",
"absences")
],
col = c("firebrick", "green3")[data_xlsx_G3_cathegorical$G3], pch = 19)
Vemos que al disponer de un número de variables elevado, no se aprecia apenas el gráfico, por lo que no podemos deducir a primera vista alguna pareja de variables que separe adecuadamente las calificaciones. Igualmente, no parece que la haya.
5.2 Normalidad univariante y multivariante
5.2.1 Normalidad univariante
A continuación se realiza una exploración gráfica de la normalidad de las distribuciones univariantes de nuestros predictores representando los histogramas y los gráficos qqplots.
Histogramas univariantes
# Histogram representation of each variable for each species
par(mfcol = c(2, 3))
for (k in 1:18) {
j0 <- names(train_data)[k]
x0 <- seq(min(train_data[, k]), max(train_data[, k]), le = 50)
for (i in 1:2) {
i0 <- levels(train_data$G3)[i]
x <- train_data[train_data$G3 == i0, j0]
hist(x, proba = T, col = grey(0.8), main = paste("G3", i0), xlab = j0)
lines(x0, dnorm(x0, mean(x), sd(x)), col = "red", lwd = 2)
}
}
Gráficos qqplots
# Representation of normal quantiles of each variable for each species
par(mfcol = c(2, 3))
for (k in 1:18) {
j0 <- names(train_data)[k]
x0 <- seq(min(train_data[, k]), max(train_data[, k]), le = 50)
for (i in 1:2) {
i0 <- levels(train_data$G3)[i]
x <- train_data[train_data$G3 == i0, j0]
qqnorm(x, main = paste("G3", i0, j0), pch = 19, col = i + 1)
qqline(x)
}
}
Este análisis exploratorio puede darnos una idea de la posible distribución normal de las variables univariadas, pero siempre es mejor hacer los respectivos test de normalidad.
Test de normalidad univariantes (Shapiro-Wilks)
Se contrasta la hipótesis nula de que los datos siguen una distribución normal univariante. Se rechaza esta hipótesis si el p-valor dado por el test de Shapiro-Wilk es menor que 0.05. En otro caso no se rechaza el supuesto de normalidad de los datos.
data_tidy <- melt(train_data, value.name = "valor")
aggregate(valor ~ G3 + variable, data = data_tidy,
FUN = function(x){shapiro.test(x)$p.value})Se encuentran evidencias de falta de normalidad univariante (p-valor < 0.05). A pesar de ello, se procede como en ACP y AF, es decir, tiramos para adelante con el fin de mostrar cómo se realizaría el análisis discriminante, pero habría que tener cuidado a la hora de obtener los resultados.
5.2.2 Normalidad multivariante
Test de normalidad multivariante (Mardia, Henze-Zirkler y Royston)
El paquete MVN contiene funciones que permiten realizar los tres test que se utilizan habitualmente para contrastar la normalidad multivariante.
Esta normalidad multivariante puede verse afectada por la presencia de outliers multivariantes. En este paquete también encontramos funciones para el análisis de outliers multivariantes. Se haría de la siguiente forma:
\(outliers <- mvn(data = train\_data[, -ncol(train\_data)], mvnTest = "hz", multivariateOutlierMethod = "quan")\)
Si la ejecutamos obtenemos el error de que la matriz de correlaciones no tiene inversa ya que tiene determinante igual a cero, lo cual no es cierto como se puede ver a continuación.
cor_m <- cor(train_data[, -ncol(train_data)])
det(cor_m)## [1] 0.2048836De igual forma, los dos test realizados a continuación encuentran evidencias al 5% de significación de falta de normalidad multivariante.
# Royston multivariate normality test
royston_test <- mvn(data = train_data[,-ncol(train_data)], mvnTest = "royston", multivariatePlot = "qq")
royston_test$multivariateNormality# Henze-Zirkler multivariate normality test
hz_test <- mvn(data = train_data[,-ncol(train_data)], mvnTest = "hz")
hz_test$multivariateNormalityCabe destacar que el discriminante lineal es muy sensible a la falta de normalidad multivariante y habría que leer sus resultados con cuidado, pero el cuadrático es bastante robusto frente a esta ausencia.
5.3 Homogeneidad de la varianza
Para el análisis de la homogeneidad de la varianza:
Cuando hay un solo predictor el test más recomendable es el test de Barttlet.
Cuando se emplean múltiples predictores, se tiene que contrastar que la matriz de covarianzas es constante en todos los grupos. En este caso es también recomendable comprobar la homogeneidad de la varianza para cada predictor a nivel individual. El test más recomendable es el test de Box M, que es una extensión del de Barttlet para escenarios multivariantes. Hay que tener en cuenta que es muy sensible a que los datos efectivamente se distribuyan según una normal multivariante. Por este motivo se recomienda utilizar una significación (p-value) < 0.001 para rechazar la hipótesis nula.
La hipóstesis nula a contrastar es la de igualdad de matrices de covarianzas en todos los grupos.
boxM(data = train_data[, -19], grouping = train_data[, 19])##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: train_data[, -19]
## Chi-Sq (approx.) = 174.12, df = 171, p-value = 0.4192En este caso no se rechaza la hipótesis nula ya que p-valor > 0.001 y por tanto se asume la homogeneidad de varianzas. Es importante recordar que para que esta conclusión sea fiable debe darse el supuesto de normalidad multivariante (la cual vimos que no se daba en nuestro caso, pero aún así tiramos para adelante).
5.4 Función discriminante
Aunque no se satisfagan los supuestos de normalidad multivariante, vamos a ajustar un modelo de clasificación discriminante lineal (tiramos para adelante con el fin de mostrar cómo se realizaría el análisis discriminante lineal, aunque hay que tenerlo presente ante la posiblidad de obtener resultados inesperados).
La función lda() del paquete MASS ajusta este modelo.
lda_model <- lda(formula = G3 ~ sex + age + famsize + Mjob + reason + traveltime +
studytime + famsup + paid + activities + romantic +
famrel + freetime + goout + Dalc + Walc + health +
absences,
data = train_data)
lda_model## Call:
## lda(G3 ~ sex + age + famsize + Mjob + reason + traveltime + studytime +
## famsup + paid + activities + romantic + famrel + freetime +
## goout + Dalc + Walc + health + absences, data = train_data)
##
## Prior probabilities of groups:
## aprobado suspenso
## 0.6708861 0.3291139
##
## Group means:
## sex age famsize Mjob reason traveltime studytime
## aprobado 0.4952830 16.55660 0.7216981 2.469297 1.301887 1.377249 2.073822
## suspenso 0.4615385 17.04503 0.7211538 2.721154 1.288462 1.405635 2.070545
## famsup paid activities romantic famrel freetime goout
## aprobado 0.3962264 0.5141509 0.4811321 0.6981132 4.157258 3.349483 3.009434
## suspenso 0.3557692 0.6442308 0.5000000 0.5769231 4.006931 3.395100 3.432692
## Dalc Walc health absences
## aprobado 1.336595 2.330189 3.551887 3.907037
## suspenso 1.432896 2.480769 3.740385 4.229504
##
## Coefficients of linear discriminants:
## LD1
## sex -0.158424951
## age 0.387186233
## famsize -0.015049772
## Mjob 0.162799150
## reason -0.048869654
## traveltime 0.069702592
## studytime -0.882737708
## famsup -0.711541738
## paid 0.850381465
## activities 0.003553655
## romantic -0.529534163
## famrel -0.543551569
## freetime -0.065047956
## goout 0.517370976
## Dalc 0.370985126
## Walc -0.114904637
## health 0.193706203
## absences -0.006099256La salida de este objeto, muestra las probabilidades a priori de cada grupo, en este caso 0.6708861 para aprobado y 0.3291139 para suspenso, las medias de cada regresor por grupo y los coeficientes del modelo de clasificación discriminante lineal.
Una vez construido el clasificador, se pueden clasificar nuevos datos en función de sus medidas sin más que llamar a la función predict(). Por ejemplo, vamos a clasificar todas las observaciones del dataset test.
new_observations <- test_data
predict(object = lda_model, newdata = new_observations)## $class
## [1] aprobado aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [9] suspenso aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [17] aprobado aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [25] aprobado aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [33] aprobado aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [41] aprobado aprobado aprobado aprobado aprobado suspenso aprobado aprobado
## [49] aprobado aprobado aprobado aprobado aprobado suspenso aprobado aprobado
## [57] aprobado aprobado aprobado aprobado aprobado aprobado suspenso aprobado
## [65] aprobado suspenso aprobado aprobado aprobado aprobado aprobado aprobado
## [73] suspenso aprobado aprobado aprobado aprobado aprobado aprobado
## Levels: aprobado suspenso
##
## $posterior
## aprobado suspenso
## 1 0.5287203 0.47127972
## 12 0.8322448 0.16775522
## 21 0.9583155 0.04168449
## 27 0.8581690 0.14183098
## 29 0.7140187 0.28598133
## 30 0.6233143 0.37668572
## 32 0.8305574 0.16944257
## 37 0.8448085 0.15519145
## 41 0.3992119 0.60078808
## 50 0.6752142 0.32478580
## 51 0.7367466 0.26325341
## 55 0.9487876 0.05121244
## 56 0.8248841 0.17511593
## 57 0.9180837 0.08191632
## 59 0.7584683 0.24153170
## 77 0.8855864 0.11441364
## 86 0.7341848 0.26581517
## 92 0.8983315 0.10166847
## 96 0.8713945 0.12860552
## 100 0.6920912 0.30790881
## 102 0.6052546 0.39474545
## 104 0.7184537 0.28154634
## 106 0.7932028 0.20679716
## 107 0.9241282 0.07587181
## 113 0.6925321 0.30746786
## 114 0.9048116 0.09518845
## 115 0.8839749 0.11602513
## 117 0.7371635 0.26283651
## 124 0.5809519 0.41904815
## 131 0.6136074 0.38639265
## 137 0.5518007 0.44819929
## 139 0.6277175 0.37228252
## 141 0.8569860 0.14301402
## 143 0.8899853 0.11001465
## 152 0.5356066 0.46439336
## 156 0.9037751 0.09622488
## 160 0.5052616 0.49473838
## 166 0.8369355 0.16306454
## 176 0.7931131 0.20688693
## 185 0.8789448 0.12105515
## 197 0.7963019 0.20369812
## 199 0.5084081 0.49159186
## 204 0.7727307 0.22726929
## 210 0.7688319 0.23116814
## 211 0.6156840 0.38431599
## 221 0.2935173 0.70648271
## 223 0.9015213 0.09847869
## 237 0.8348371 0.16516295
## 257 0.7818562 0.21814376
## 260 0.8780754 0.12192462
## 263 0.8102649 0.18973507
## 265 0.5543106 0.44568938
## 275 0.6878363 0.31216375
## 276 0.4653202 0.53467976
## 280 0.8891487 0.11085128
## 286 0.8087309 0.19126909
## 288 0.7580155 0.24198453
## 289 0.8025402 0.19745977
## 294 0.7617490 0.23825104
## 295 0.8384137 0.16158628
## 296 0.7966999 0.20330005
## 297 0.5492399 0.45076007
## 298 0.3752390 0.62476097
## 300 0.8092060 0.19079400
## 302 0.8499428 0.15005725
## 310 0.4715284 0.52847161
## 316 0.5394510 0.46054901
## 329 0.8247123 0.17528768
## 339 0.7196361 0.28036391
## 340 0.6953299 0.30467010
## 343 0.5568235 0.44317647
## 344 0.5309468 0.46905321
## 354 0.1935134 0.80648664
## 358 0.5618232 0.43817679
## 365 0.6670069 0.33299308
## 367 0.7672386 0.23276139
## 379 0.7346458 0.26535417
## 388 0.7515503 0.24844974
## 389 0.7183630 0.28163704
##
## $x
## LD1
## 1 0.886774586
## 12 -0.982390749
## 21 -2.910417203
## 27 -1.232026790
## 29 -0.119059037
## 30 0.398135389
## 32 -0.967255462
## 37 -1.099108605
## 41 1.545311460
## 50 0.111181644
## 51 -0.262576290
## 55 -2.639027832
## 56 -0.917228217
## 57 -2.007077792
## 59 -0.407387406
## 77 -1.541663112
## 86 -0.246020608
## 92 -1.708124085
## 96 -1.374331598
## 100 0.013046932
## 102 0.493983864
## 104 -0.146496000
## 106 -0.658904361
## 107 -2.111706963
## 113 0.010444219
## 114 -1.799967316
## 115 -1.521787377
## 117 -0.265280286
## 124 0.620629736
## 131 0.449860110
## 137 0.769916606
## 139 0.374500688
## 141 -1.219848124
## 143 -1.597189067
## 152 0.851996350
## 156 -1.784910111
## 160 1.004914006
## 166 -1.025115103
## 176 -0.658216330
## 185 -1.461251985
## 197 -0.682792002
## 199 0.989086100
## 204 -0.507338589
## 210 -0.479591758
## 211 0.438836347
## 221 2.135769816
## 223 -1.752660679
## 237 -1.005881887
## 257 -0.573627127
## 260 -1.451009218
## 263 -0.793931666
## 265 0.757149679
## 275 0.038056697
## 276 1.206073713
## 280 -1.586481128
## 286 -0.781424235
## 288 -0.404281403
## 289 -0.731712057
## 294 -0.430009130
## 295 -1.038784224
## 296 -0.685879859
## 297 0.782928478
## 298 1.672372407
## 300 -0.785289623
## 302 -1.149026703
## 310 1.174725459
## 316 0.832552201
## 329 -0.915733887
## 339 -0.153855209
## 340 -0.006118322
## 343 0.744353355
## 344 0.875536714
## 354 2.826017934
## 358 0.718849065
## 365 0.157935977
## 367 -0.468347528
## 379 -0.248992315
## 388 -0.360359713
## 389 -0.145932280Por ejemplo, según la función discriminante, la probabilidad a posteriori de que el segundo registro del conjunto test apruebe es del 89.6% mientras la de que suspenda es tan solo del 10.4%. Por tanto este dato sería clasificado en la calificación aprobado. Se razona del mismo modo con el resto de elementos del conjunto test.
5.5 Validación del modelo
La función confusionmatrix() del paquete biotools realiza una validación cruzada del modelo de clasificación.
pred <- predict(object = lda_model, newdata = test_data)
confusionmatrix(test_data$G3, pred$class)## new aprobado new suspenso
## aprobado 50 3
## suspenso 23 3# Classification error percentage
test_error <- mean(test_data$G3 != pred$class) * 100
paste("test_error=", test_error, "%")## [1] "test_error= 32.9113924050633 %"5.6 Visualización de las clasificaciones
Desde el punto de vista geométrico, el análisis discriminante lineal separa el espacio mediante una recta. En este sentido, la función partimat() del paquete klaR permite representar los límites de clasificación de un modelo discriminante lineal o cuadrático para cada par de predictores. Cada color representa una región de clasificación acorde al modelo, se muestra el centroide de cada región y el valor real de las observaciones.
partimat(G3 ~ sex + age + famsize + Mjob + reason + traveltime + studytime +
famsup + paid + activities + romantic + famrel + freetime +
goout + Dalc + Walc + health + absences,
data = test_data, method = "lda", prec = 200,
image.colors = c("green", "orange"),
col.mean = "yellow",nplots.vert =1, nplots.hor=3)
6 Análisis Discriminante Cuadrático (ADC)
Al igual que para el Análisis Discriminante Lineal, para realizar el Análisis Discriminante Cuadrático se empieza con la exploración gráfica de los datos y las comprobaciones sobre normalidad univariante y multivariante y homogeneidad de las varianzas, que ya se han realizado previamente.
Vimos que no se verifica el supuesto de normalidad multivariante. Aun así, se procede a ajustar un modelo discriminante cuadrático porque es robusto frente a la falta de este supuesto, aunque hay que tenerlo presente ante la posiblidad de obtener resultados inesperados.
6.1 Función discriminante
La función qda() del paquete MASS realiza la clasificación.
qda_model <- qda(formula = G3 ~ sex + age + famsize + Mjob + reason + traveltime +
studytime + famsup + paid + activities + romantic +
famrel + freetime + goout + Dalc + Walc + health +
absences,
data = train_data)
qda_model## Call:
## qda(G3 ~ sex + age + famsize + Mjob + reason + traveltime + studytime +
## famsup + paid + activities + romantic + famrel + freetime +
## goout + Dalc + Walc + health + absences, data = train_data)
##
## Prior probabilities of groups:
## aprobado suspenso
## 0.6708861 0.3291139
##
## Group means:
## sex age famsize Mjob reason traveltime studytime
## aprobado 0.4952830 16.55660 0.7216981 2.469297 1.301887 1.377249 2.073822
## suspenso 0.4615385 17.04503 0.7211538 2.721154 1.288462 1.405635 2.070545
## famsup paid activities romantic famrel freetime goout
## aprobado 0.3962264 0.5141509 0.4811321 0.6981132 4.157258 3.349483 3.009434
## suspenso 0.3557692 0.6442308 0.5000000 0.5769231 4.006931 3.395100 3.432692
## Dalc Walc health absences
## aprobado 1.336595 2.330189 3.551887 3.907037
## suspenso 1.432896 2.480769 3.740385 4.229504La salida de este objeto, muestra las probabilidades a priori de cada grupo, en este caso 0.6708861 para aprobado y 0.3291139 para suspenso, las medias de cada regresor por grupo.
Una vez construido el clasificador, se pueden clasificar nuevos datos en función de sus medidas sin más que llamar a la función predict(). Por ejemplo, vamos a clasificar todas las observaciones del dataset test.
new_observations <- test_data
predict(object = qda_model, newdata = new_observations)## $class
## [1] aprobado aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [9] suspenso suspenso aprobado aprobado aprobado aprobado aprobado aprobado
## [17] aprobado suspenso aprobado aprobado aprobado suspenso aprobado aprobado
## [25] aprobado aprobado aprobado aprobado suspenso suspenso aprobado aprobado
## [33] aprobado aprobado aprobado aprobado suspenso aprobado aprobado suspenso
## [41] aprobado aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [49] aprobado aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [57] aprobado aprobado aprobado aprobado aprobado aprobado aprobado aprobado
## [65] aprobado aprobado aprobado aprobado aprobado suspenso aprobado suspenso
## [73] suspenso aprobado aprobado aprobado aprobado aprobado aprobado
## Levels: aprobado suspenso
##
## $posterior
## aprobado suspenso
## 1 0.9420895 5.791050e-02
## 12 0.9132794 8.672064e-02
## 21 0.9983882 1.611764e-03
## 27 0.9996942 3.058088e-04
## 29 0.9519198 4.808024e-02
## 30 0.8740865 1.259135e-01
## 32 0.9725168 2.748321e-02
## 37 0.9933374 6.662618e-03
## 41 0.1599467 8.400533e-01
## 50 0.3660455 6.339545e-01
## 51 0.9921764 7.823598e-03
## 55 0.9978306 2.169449e-03
## 56 0.9628626 3.713736e-02
## 57 0.9866293 1.337075e-02
## 59 0.9659756 3.402435e-02
## 77 0.7117161 2.882839e-01
## 86 0.6263690 3.736310e-01
## 92 0.4917608 5.082392e-01
## 96 0.7509006 2.490994e-01
## 100 0.8578510 1.421490e-01
## 102 0.6992017 3.007983e-01
## 104 0.2780915 7.219085e-01
## 106 0.9881970 1.180300e-02
## 107 0.9652249 3.477514e-02
## 113 0.5679399 4.320601e-01
## 114 0.9961207 3.879346e-03
## 115 0.9992368 7.632085e-04
## 117 0.9881334 1.186660e-02
## 124 0.1060614 8.939386e-01
## 131 0.1369429 8.630571e-01
## 137 0.9139531 8.604690e-02
## 139 0.5283120 4.716880e-01
## 141 0.7562291 2.437709e-01
## 143 0.9327727 6.722733e-02
## 152 0.6146317 3.853683e-01
## 156 0.9999709 2.906058e-05
## 160 0.4715576 5.284424e-01
## 166 0.8308611 1.691389e-01
## 176 0.6283201 3.716799e-01
## 185 0.3972366 6.027634e-01
## 197 0.9960218 3.978206e-03
## 199 0.8719620 1.280380e-01
## 204 0.9873704 1.262960e-02
## 210 0.7181323 2.818677e-01
## 211 0.5816051 4.183949e-01
## 221 0.5186015 4.813985e-01
## 223 0.9791096 2.089043e-02
## 237 0.8131715 1.868285e-01
## 257 0.7129042 2.870958e-01
## 260 0.9856287 1.437128e-02
## 263 0.9499050 5.009498e-02
## 265 0.7846896 2.153104e-01
## 275 0.7546396 2.453604e-01
## 276 0.9453512 5.464877e-02
## 280 0.9917348 8.265192e-03
## 286 0.9869410 1.305902e-02
## 288 0.5607798 4.392202e-01
## 289 0.9415090 5.849101e-02
## 294 0.6537919 3.462081e-01
## 295 0.8604681 1.395319e-01
## 296 0.6809333 3.190667e-01
## 297 0.9434945 5.650552e-02
## 298 0.7037053 2.962947e-01
## 300 0.8192494 1.807506e-01
## 302 0.8940954 1.059046e-01
## 310 0.8961285 1.038715e-01
## 316 0.5901816 4.098184e-01
## 329 0.9292400 7.076001e-02
## 339 0.9845709 1.542907e-02
## 340 0.4540543 5.459457e-01
## 343 0.8105703 1.894297e-01
## 344 0.4402601 5.597399e-01
## 354 0.1499105 8.500895e-01
## 358 0.8961398 1.038602e-01
## 365 0.8286849 1.713151e-01
## 367 0.9680180 3.198197e-02
## 379 0.6149508 3.850492e-01
## 388 0.9583831 4.161691e-02
## 389 0.9928296 7.170400e-03Por ejemplo, según la función discriminante, la probabilidad a posteriori de que el segundo registro del conjunto test apruebe es del 99.9% mientras la de que suspenda es tan solo del 0.1%. Por tanto este dato sería clasificado en la calificación aprobado. Se razona del mismo modo con el resto de elementos del conjunto test.
6.2 Validación del modelo
La función confusionmatrix() del paquete biotools realiza una validación cruzada del modelo de clasificación.
pred <- predict(object = qda_model, newdata = test_data)
confusionmatrix(test_data$G3, pred$class)## new aprobado new suspenso
## aprobado 47 6
## suspenso 21 5# Classification error percentage
test_error <- mean(test_data$G3 != pred$class) * 100
paste("test_error=", test_error, "%")## [1] "test_error= 34.1772151898734 %"6.3 Visualización de las clasificaciones
La función partimat() del paquete klaR permite representar los límites de clasificación de un modelo discriminante lineal o cuadrático para cada par de predictores. Cada color representa una región de clasificación acorde al modelo, se muestra el centroide de cada región y el valor real de las observaciones.
partimat(G3 ~ sex + age + famsize + Mjob + reason + traveltime + studytime +
famsup + paid + activities + romantic + famrel + freetime +
goout + Dalc + Walc + health + absences,
data = test_data, method = "qda", prec = 200,
image.colors = c("darkgoldenrod1", "skyblue2"),
col.mean = "firebrick", nplots.vert =1, nplots.hor=3)
LDA vs QDA
El clasificador más adecuado depende de las implicaciones que tenga asumir que todos los grupos comparten una matriz de covarianza común ya que este supuesto puede producir un sesgo en las clasificaciones o producir varianzas altas.
- LDA produce límites de decisión lineales, lo que se traduce en menor flexibilidad y por lo tanto menor problema de varianza.
- QDA produce límites cuadráticos y por lo tanto curvos, lo que aporta mayor flexibilidad permitiendo ajustarse mejor a los datos, menor sesgo pero mayor riesgo de varianza.
- En términos generales, LDA tiende a conseguir mejores clasificaciones que QDA cuando hay pocas observaciones con las que entrenar al modelo, escenario en el que evitar la varianza es crucial.
- Si se dispone de una gran cantidad de observaciones de entrenamiento (como es nuestro caso) o si no es asumible que existe una matriz de covarianza común entre clases, QDA es más adecuado.
7 Análisis cluster (AC)
Vamos a realizar un Análisis Clúster que confirme que la agrupación de la variable respuesta (variable G3 categórica) utilizada en los modelos de clasificación lineal y cuadrático es adecuada.
7.1 Algunas distancias para Análisis Cluster
Para clasificar las observaciones en grupos es necesario elegir medidas de similaridad, o de distancia (disimilaridad), adecuadas que proporcionen información de como de parecidas son dos observaciones cualesquiera. De hecho esta elección influye en el tamaño y forma de los clusters. Esta elección es un paso fundamental en clustering.
Algunas medidas de distancia clásicas frecuentemente utilizadas son la distancia Euclídea o la distancia Manhattan.
Otras medidas de disimilaridad ampliamente utilizadas, por ejemplo, en el análisis de datos de expresión génica son las distancias basadas en correlaciones. Estas distancias se obtienen restando la correspondiente medida de correlación al valor 1. Entre estas medidas de distancia destacan: la distancia basada en la correlación de Pearson, en la correlación de Spearman, correlación de Kendall, etc.
Como se ha dicho anteriormente la elección de la distancia es muy importante. Casi todo el software habitual para Análisis Cluster utiliza la distancia euclídea, aunque dependiendo del tipo de datos y de las preguntas de investigación planteadas puede interesar otra medida de disimilaridad o distancia.
En R es fácil calcular y visualizar la matriz de distancias entre observaciones con las funciones get_dist() y fviz_dist(), respectivamente, incluidas en el paquete factoextra. De hecho, aunque por defecto, get_dist() calcula la distancia euclídea también admite como parámetro todas las distancias mencionadas anteriormente.
La siguiente matriz de distancias muestra en marrón aquellos estados que presentan grandes disimilitudes (distancias), frente a aquellos que parecen más cercanos en amarillo. Se utiliza el color blanco para referirse a aquellos estados con distancias no tan extremas como para ser consideradas como bajas o altas.
distance<- get_dist(data_xlsx_reduced)
fviz_dist(distance, gradient = list(low ="yellow", mid = "white", high = "brown"))
7.2 Clustering jerárquico: método de Ward
El agrupamiento jerárquico está interesado en encontrar una jerarquía basada en la cercanía o semejanza de los datos según la distancia considerada. En el caso aglomerativo, se parte de un grupo con las observaciones más cercanas. A continuación se calculan los siguientes pares más cercanos y de manera ascendente se van generando grupos. Esta construcción se puede observar de forma visual con la construcción de un dendrograma.
A continuación se ilustrará cómo los grupos están definidos por la cantidad de líneas verticales del dendrograma, y la selección del número de grupos óptimo se podrá estimar desde este mismo gráfico.
dendrogram <- hclust(dist(data_xlsx_reduced, method = 'euclidean'), method = 'ward.D')
ggdendrogram(dendrogram, rotate = FALSE, labels = FALSE, theme_dendro = TRUE) +
labs(title = "Dendrograma")
En el eje horizontal del dendrograma tenemos cada uno de los datos que componen el conjunto de entrada, mientras que en el eje vertical se representa la distancia euclídea que existe entre cada grupo a medida que éstos se van jerarquizando.
Cada línea vertical del diagrama representa un agrupamiento. Conforme se va subiendo en el dendograma termina con un solo gran grupo determinado por la línea horizontal superior. De modo que, al ir descendiendo en la jerarquía, se observa que de un solo grupo pasamos a 2, luego a 3, luego a 4, y así sucesivamente.
Una forma de determinar el número K de grupos adecuado es cortando el dendrograma a aquella altura del diagrama que mejor representa los datos de entrada.
7.3 Clustering no jerárquico: algoritmo K-means
R implementa el algoritmo K-means con la función del mismo nombre. Esta función recibe como parámetros de entrada los datos y el número de agrupamientos a realizar (parámetro centers). Para abordar el problema de la elección de los puntos semilla iniciales incorpora el parámetro nstart que prueba múltiples configuraciones iniciales e informa sobre la mejor. Por ejemplo, si nstart = 25, generará 25 configuraciones iniciales. El uso de este parámetro es recomendable.
Para este primer ejemplo la función kmeans() construye dos clusters.
k2 <- kmeans(data_xlsx_reduced, centers = 2, nstart = 25)
# Displaying all the fields of the object k2
str(k2)## List of 9
## $ cluster : int [1:395] 1 1 1 1 1 1 1 1 1 1 ...
## $ centers : num [1:2, 1:19] 0.484 0.386 16.646 16.977 0.698 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:2] "1" "2"
## .. ..$ : chr [1:19] "sex" "age" "famsize" "Mjob" ...
## $ totss : num 19315
## $ withinss : num [1:2] 13240 750
## $ tot.withinss: num 13990
## $ betweenss : num 5325
## $ size : int [1:2] 351 44
## $ iter : int 1
## $ ifault : int 0
## - attr(*, "class")= chr "kmeans"La salida que proporciona la función kmeans() es una lista de información, de la que destacan las siguientes:
- cluster: es un vector de enteros, de 1 a K (K = 2 en este caso), que indica el cluster en el que ha sido ubicado cada observación.
- centers: una matriz con los sucesivos centros de los clusters.
- totss: la suma total de cuadrados.
- withinss: vector de suma de cuadrados dentro de cada cluster (un componente por cluster).
- tot.withinss: suma total de cuadrados de los cluster, i.e. sum(withinss).
- betweens: suma de cuadrados entre grupos, i.e. totss-tot.withinss.
- size: el número de observaciones en cada cluster.
Al mostrar la variable k2 se ve como las agrupaciones dan como resultado 2 tamaños de agrupación de 44 y 351 También se ven los centros de cada grupo (medias) en todas las variables. Y por último la asignación de grupo para cada observación.
k2## K-means clustering with 2 clusters of sizes 351, 44
##
## Cluster means:
## sex age famsize Mjob reason traveltime studytime famsup
## 1 0.4843305 16.64582 0.6980057 2.458652 1.270655 1.379977 2.072080 0.3874644
## 2 0.3863636 16.97727 0.8181818 2.818182 1.295455 1.517970 2.076143 0.3863636
## paid activities romantic famrel freetime goout Dalc Walc
## 1 0.5099715 0.4985755 0.6894587 4.109019 3.353057 3.079772 1.333515 2.313390
## 2 0.7954545 0.4318182 0.4772727 4.054816 3.311291 3.340909 1.363636 2.113636
## health absences G3
## 1 3.535613 4.4491379 11.62678
## 2 3.704545 0.2860918 0.75000
##
## Clustering vector:
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1
## [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1
## [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 2 1 1 2 2 2 2 1 1 2 1 1 1 2 1 2 1
## [149] 2 1 2 1 1 2 1 1 1 1 1 1 2 1 2 1 2 1 1 1 2 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 2
## [223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [260] 2 1 1 1 1 2 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [297] 2 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2
## [334] 2 2 1 1 2 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1
## [371] 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 1 2 1 1 1 1 1
##
## Within cluster sum of squares by cluster:
## [1] 13240.3490 750.1396
## (between_SS / total_SS = 27.6 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"Una forma visual de resumir los resultados de forma elegante y con una interpretación directa es mediante el uso de la función fviz_cluster().
fviz_cluster(k2, data=data_xlsx_reduced)
Observación: Si hay más de dos dimensiones (variables), esta función realizará, en primer lugar, un análisis de componentes principales (ACP) y dibujará los puntos de acuerdo a las dos primeras componentes principales obtenidas (las que explican la mayor parte de la varianza). Es por esto que en el gráfico anterior, Dim1 y Dim2 se refieren a estas dos componentes principales.
Para usar el algoritmo K-means, el número K de clusters debe ser fijado de antemano. Es por esto que es recomendable ejecutar el mismo proceso con otros valores de K (en este ejemplo para K = 3, 4 y 5) para comparar y examinar las diferencias entre los resultados.
set.seed(123)
k3 <- kmeans(data_xlsx_reduced, centers = 3, nstart = 25)
k4 <- kmeans(data_xlsx_reduced, centers = 4, nstart = 25)
k5 <- kmeans(data_xlsx_reduced, centers = 5, nstart = 25)
# Plots to compare
p24 <- fviz_cluster(k2, geom = "point", data = data_xlsx_reduced) + ggtitle("k = 2")
p25 <- fviz_cluster(k3, geom = "point", data = data_xlsx_reduced) + ggtitle("k = 3")
p26 <- fviz_cluster(k4, geom = "point", data = data_xlsx_reduced) + ggtitle("k = 4")
p27 <- fviz_cluster(k5, geom = "point", data = data_xlsx_reduced) + ggtitle("k = 5")
grid.arrange(p24, p25, p26, p27, nrow = 2)
Aunque esta visualización nos permite deducir donde ocurren las verdaderas diferencias (o no ocurren, ya que vemos que los clusters se superponen), no nos dice cuál es el número óptimo de clusters.
7.3.1 Determinación del número óptimo de clusters
Tal y como se ha indicado anteriormente, cuando se aplica un método no jerárquico como K-means para realizar un análisis cluster, el investigador debe informar a priori del número de clusters deseado. En este sentido, este investigador estará interesado en proporcionar de partida un número óptimo de grupos a formar.
A continuación se presentan tres de los métodos más utilizados para determinar este número óptimo de grupos: método de Elbow, método de Silhouette y el stadístico Gap.
7.3.1.1 Método de Elbow
Teniendo en mente que la idea tras una división en K clusters, es obtener estas agrupaciones de modo que la varianza total intra-grupos sea mínima (total within-cluster variation o total within-cluster sum of square), se puede utilizar el siguiente algoritmo para identificar el número óptimo de clusters:
- Ejecutar un algoritmo de clustering (como K-means) para diferentes valores del K (por ejemplo K = 1,…,10).
- Para cada K se calcula la variación total intra-cluster (total within-cluster sum of square, que aquí denotamos por wss).
- Se dibuja la curva de wss de acuerdo al número de clusters K. La localización en esta curva de un “codo” es tomado como el indicador más apropiado del número de clusters.
Aunque podemos programar este algoritmo en R, el método de Elbow está implementado en la función fviz_nbclust().
set.seed(123)
fviz_nbclust(data_xlsx_reduced, kmeans, method = "wss")
7.3.1.2 Método de Silhouette
Este enfoque, silueta promedio, mide la calidad de un cluster. Es decir, determina como de adecuado es un objeto dentro de su cluster. El número óptimo de clusters según este enfoque es, de entre un rango de valores posibles para K, aquel que maximiza la silueta promedio.
Como antes, este algoritmo se puede programar en R, pero la función fviz_nbclust() también lo incluye.
set.seed(123)
fviz_nbclust(data_xlsx_reduced, kmeans, method = "silhouette")
7.3.1.3 Método estadístico de brecha (GAP)
Este método compara la variación total intracluster para diferentes valores de K. Utiliza simulación Montecarlo en su algoritmo.
La función clusGap() proporciona el estadístico GAP y su error estándar para una salida. Con la función fviz_gap_stat() se obtiene una representación gráfica que sugiere un número de clusters.
set.seed(123)
gap_stat <- clusGap(data_xlsx_reduced, FUN = kmeans, nstart = 25, K.max = 10, B = 50)
fviz_gap_stat(gap_stat)
7.3.2 Análisis de los resultados
Tras el análisis pormenorizado del número óptimo de clusters, realizado en la sección anterior, y las pruebas anteriores con diferentes números de clusters, concluímos que no parece adecuado realizar un Análisis Cluster. Esto puede deberse quizá a la independencia de las variables que observamos a la hora de realizar ACP o AF o la forma en la que se han codificado las variables categóricas (podría probarse a hacer este análisis cluster usando medidas de distancia para variables categóricas, como puede ser la distancia de Manhattan). Se observa que a mayor número de clusters, peor (ya que los clusters se superponen)
Aun así, si tuviesemos que realizarlo, K = 2 parece ser el número de grupos más adecuado. Esto confirma que la agrupación de la variable respuesta utilizada en los modelos de clasificación lineal y cuadrático es adecuada. Por lo tanto, el categorizar la calificación de los alumnos como aprobado o suspenso, fue una buena opción tomada a la hora de realizar los modelos de clasificación (por ejemplo, podríamos haber pensado también en categorizar la calificación de los alumnos como suspenso, aprobado, notable y sobresaliente).
set.seed(123)
final <- kmeans(data_xlsx_reduced, centers = 2, nstart = 25)
print(final)## K-means clustering with 2 clusters of sizes 351, 44
##
## Cluster means:
## sex age famsize Mjob reason traveltime studytime famsup
## 1 0.4843305 16.64582 0.6980057 2.458652 1.270655 1.379977 2.072080 0.3874644
## 2 0.3863636 16.97727 0.8181818 2.818182 1.295455 1.517970 2.076143 0.3863636
## paid activities romantic famrel freetime goout Dalc Walc
## 1 0.5099715 0.4985755 0.6894587 4.109019 3.353057 3.079772 1.333515 2.313390
## 2 0.7954545 0.4318182 0.4772727 4.054816 3.311291 3.340909 1.363636 2.113636
## health absences G3
## 1 3.535613 4.4491379 11.62678
## 2 3.704545 0.2860918 0.75000
##
## Clustering vector:
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1
## [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1
## [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 2 1 1 2 2 2 2 1 1 2 1 1 1 2 1 2 1
## [149] 2 1 2 1 1 2 1 1 1 1 1 1 2 1 2 1 2 1 1 1 2 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 2
## [223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [260] 2 1 1 1 1 2 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [297] 2 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2
## [334] 2 2 1 1 2 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1
## [371] 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 1 2 1 1 1 1 1
##
## Within cluster sum of squares by cluster:
## [1] 13240.3490 750.1396
## (between_SS / total_SS = 27.6 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"# Final output
fviz_cluster(final, data = data_xlsx_reduced)